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ABSTRACT 


A  constitutive  relationship  was  used  to  model  the  cyclic  material  response  of  damping 
test  samples  in  separate  bending  and  torsion  configurations.  This  was  done  in  order  to  better 
understand  variations  in  reported  values  of  damping  for  materials  possessing  strain  dependent 
characteristics.  The  constitutive  equations  are  based  on  a  model  of  shape  memory  alloy  stress- 
strain  behavior  and  have  been  adapted  especially  for  the  study  of  nonlinear  hysteresis  and  the 
problem  of  strain  dependent  damping,  ^perimental  measurements  and  analytical  material 
response  analyses  of  separate  bending  and  torsion  test  samples  indicated  that  when  the  damping 
of  a  single  nonlinear  material  is  plotted  against  the  one-dimensional  local  strain  of  the  sample, 
results  are  produced  which  are  difGcult  to  compare.  However,  when  the  same  results  are  plotted 
against  an  invariant  measure  of  three-dimensional  distortion  the  means  by  which  one  may 
compare  the  data  is  more  straightforward.  Also,  the  approach  allows  for  a  quantitative 
comparison  of  the  damping  at  a  material  point  to  the  overall  damping.  The  method  can  be 
applied  to  any  homogeneous  isotropic  nonlinear  damping  material. 
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0115)  Program  Element  62254N  Task  Area  RS3454,  Work  Unit  1-2812-949.  This  report 
satisfies  milestone  53SR1/4. 


INTRODUCTION 

In  addition  to  add-on  damping  techniques  currently  being  used  in  the  Navy,  material 
damping  is  being  investigated  as  a  potential  means  of  further  reducing  machine  vibration,  noise, 
and  sound  emission  in  seafaring  vessels.  Ideally  a  high  damping  structural  material  provides  a 
sufGcient  amount  of  both  stiffness  and  damping  so  as  to  be  used  as  a  sole  machine  part  or 
vibrating  element  without  added  treatments.  Such  materials  are  most  useful  for  oscillating  parts 
or  elements  that  cannot  be  damped  by  conventional  external  treatments.  Also  these  materials 
can  be  useful  in  situations  where  heat  or  other  environmental  Actors  (e.g.  moisture,  corrosion) 
have  to  be  considered.  This  approach  is  also  useful  in  damping  longitudinal  vibrations  which 
cannot  be  effectively  controlled  by  external  treatments. 
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Because  high  stiffness  and  strength  are  required  in  many  important  applications,  me^ls 
which  possess  a  large  inherent  damping  capacity  have  been  extensively  sought  [1,2].  Some 
specific  applications  include  gears  and  gear  webs,  pump  casings,  engine  parts,  propellers,  and 
others  (see  [1]).  High  damping  metals  are  also  used  as  plug  inserts  and  cladding,  and  such 
applications  can  provide  a  reduction  of  resonant  amplification  factors  as  well  as  the  attenuation 
of  ringing  [1]. 

Generally  "high  damping"  in  metal  is  a  measured  peak  loss  factor  or  phase  lag  with  a 
value  of  gr  higher.  Indeed  many  alloy  compositions  have  been  studied  and  found  to  possess 
such  levels  of  damping  (e.g.  see  [3-8]).  Mechanisms  which  give  rise  to  damping  in  metals 
include:  movement  of  point  delects,  dislocations  or  domain  walls.  These  effects  give  rise  to 
macroscopic  hysteresis  and  thus  damping.  Ihe  damping  capacity  of  high  damping  metals  is 
strain  dependent  because  the  primary  damping  mechanisms  function  over  a  finite  strain  range. 
Such  effects  give  rise  to  a  well  defined  peak  in  the  plot  of  measured  damping  vs.  specimen  strain 
amplitude.  Examples  of  magnetostrictive  metallic  materials  exhibiting  strain-dependent 
fiamping  are  given  in  Fig.  1.  This  type  of  response  is  termed  "nonlinear"  because  the  measured 
damping  capacity  varies  with  specimen  strain  amplitude. 

A  generalized  stre^-strain  diagram  corresponding  to  a  relatively  large  category  of 
nonlinear  damping  mechanisms  is  illustrated  in  Fig.  2.  Note  from  this  figure  that  the  damping 
mechanism  is  activated  near  a  critical  stress  and  becomes  saturated  at  a  strain  of  This  type 

of  hysteretic  response  can  be  as.sociated  with  a  number  of  nonlin:^r  anelastic  damping 
mechanisms.  For  example,  in  dislocation  breakaway  a  minimum  stir^ss  is  required  tc  force 
dislocations  over  nearby  pinning  points  during  loading.  Upon  unloading  the  elastic  strain  energy 
stored  in  the  lattice  of  the  material  may  be  sufficient  to  move  the  dislocations  back  to  their 
original  positions.  The  net  effect  of  this  process  is  an  elastic  response  with  internal  friction. 

Other  nonlinear  anelastic  damping  mechanisms  include  the  movement  of  mobile  domain 
boundaries.  Ferromagnetic  domain  walls,  twin  boundaries,  antiferromagnetic  domain  walls,  and 
phase  domain  walls  fall  into  this  category.  Usually,  a  finite  amount  of  stress  is  required  to 
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initiate  this  type  of  mechanism,  i.e.  a  specific  amount  of  stress  is  required  to  overcome  an  energy 
barrier  so  that  the  boundaries  may  move.  Release  of  the  applied  stress  subsequently  causes  the 
mechanism  to  act  in  reve*^  because  the  elastic  strain  energy  stored  in  the  material  is  sufficient 
to  move  the  boundaries  back  across  the  energy  barrier.  Again  the  net  effect  is  an  elastic 
response  with  energy  dissipation  due  to  internal  friction.  Also,  these  mechanisms  may  become 
saturated  at  a  limiting  value  of  strain.  For  example,  the  magnetic  domains  in  high  damping 
ferromagnetic  alloys  are  ananged  in  a  randomly  criented  pattern  when  the  material  is  unstressed. 
However,  upon  application  of  a  uniaxial  stress  the  domains  change  their  orientation  and  tend  to 
align  themselves  in  the  direction  of  loading  as  the  stress  is  increased.  Once  these  domains 
become  fully  aligned  any  further  stress  cannot  cause  relative  motion  of  the  domains  and  the 
mechanism  is  said  to  be  saturated.  The  amount  of  energy  that  can  be  dissipated  by  this  type  of 
damping  mechanism  is  therefore  limited  to  a  fixed  value  for  cyclic  strain  amplitudes  greater  than 
the  limiting  value  of  strain  corresponding  to  saturation  (e^). 

The  data  obtained  for  a  single  strain  dependent  material  in  varied  test  configurations  is 
often  difficult  to  compare  because  of  the  inherent  strain  distributions  that  arise  from  the  loading. 
Data  from  separate  bending  and  torsion  tests  [3]  given  in  Fig.  3  shows  this  effect;  indeed  the 
results  indicate  that  the  torsional  tests  produce  significantly  higher  values  of  damping  for 
common  levels  of  peak  sample  strain.  Also,  the  damping  rises  more  steeply  with  peak  sample 
strain  in  the  torsion  test.  However  it  is  important  to  note  that  the  strains  on  the  abscissa  are  shear 
strains  in  the  case  of  torsional  data  and  axial  strains  in  the  case  of  bending  data  and  these 
separate  strains  are  not  equivalent;  indeed  axial  strains  give  a  measure  of  length  change  while 
shear  strains  refer  to  the  distortion  of  right  angles  This  is  an  important  aspect  of  the  problem 
which  will  be  discussed  in  the  analysis  section  of  this  paper. 

Strain  dependent  materials  are,  at  best,  difficult  to  model  analytically  because  of  their 
nonlinear  characteristics.  Early  work  in  this  area  concentrated  on  evaluating  the  damping  of 
members  by  combining  material  energy  absorbing  properties  with  geometric  and  stress 
distribution  factors  [9,10].  Another  approach  is  to  use  a  constitutive  law  which  describes 


DTRC-SME-91/34 


3 


nonlinear  material  behavior  and  hysteresis  at  a  point,  and  this  approach  will  be  used  here.  Many 
such  laws  exist  (e.g.  see  [11,12]),  but  these  are  usually  specific  to  postyielding  viscoplastic 
behavior  and  large  strain  levels.  In  this  paper  a  proposed  constitutive  law  [13,14]  for  the  stress- 
strain  behavior  of  shape  memory  alloys  is  adapted  to  the  case  of  nonlinear  damping.  The 
equations  of  this  law  were  applied  to  the  cases  of  simple  uniaxial  tension-compression  and  shear 
loading  of  materials.  Solid  geometries  of  beam  and  shaft  test  samples  in  bending  and  torsion 
were  also  considered.  The  strain  dependent  nature  of  each  test  configuration  was  computed,  and 
because  this  behavior  was  of  primary  interest,  temperature  and  fi^equency  effects  were  not 
considered. 

ANALYSIS 

A  three-dimensional  constitutive  law  of  hysteretic  material  behavior  was  employed  so 
that  a  useful  study  of  strain  dependent  damping  could  be  made.  This  law  is  based  on  the  three 
dimensional  generalization  [14]  of  a  one-dimensional  model  of  shape  memory  alloy  (SMA) 
stress-strain  behavior  [13],  where  the  extension  from  one  to  three  dimensions  follows  a  method 
originally  developed  by  Prager  [15]  (for  a  detailed  development  regarding  this  extension  method 
see  [16,17,14]).  This  choice  of  modeling  schemes  was  pursued  because  the  hysteretic  response 
of  superelastic  SMA's  is  very  similar  in  character  to  that  of  high  damping  metals  (see  Fig.  2), 
except  that  the  stress  and  strain  levels  are  different  by  many  orders  of  magnitude.  This  does  not 
prevent  the  use  of  the  constitutive  law,  however,  as  long  as  the  material  properties  of  the  law  can 
be  scaled  to  accommodate  the  lower  stress  and  strain  levels  associated  with  the  dissipative 
mechanisms  of  the  damping  material. 

The  constitutive  law  is  for  homogeneous  and  isotropic  material  behavior  and  is  based 
upon  a  separation  of  strain  and  strain  rate  into  elastic  and  inelastic  components: 


el  m 
y  ij  ij 


(la) 


.  .el  .10 
y  ij  y 


(lb) 


DTRC-SME-91/34 


Here  an  ove±ead  dot  represents  differentiation  with  respect  to  time.  Thus  Ejj  and  e^j  are  the 
three-dimensional  tensors  of  strain  and  strain  rate,  and  the  superscripts  "el"  and  "in"  designate 
the  respective  elastic  and  inelastic  components  of  each.  The  elastic  component  follows  directly 
from  the  theory  of  isotropic  elasticity  [17]: 


cl  1+v  V 

“  E  °ij  ■  E  ‘’I* 


(2) 


where  Ojj  is  the  stress  tensor,  6jj  is  the  Kronecker  delta^  and  where  E  and  v  are  the  elastic 
material  constants. 

The  basic  equations  for  the  evolution  of  inelastic  strain  were  taken  from  the  previously 
cited  model  of  shape  memory  alloy  behavior.  In  this  model  the  growth  of  inelastic  strain  is  a 
function  of  a  backstress  tensor  Pjj,  which  is  a  variable  that  accounts  for  internal  stress  fields  in 

the  material,  as  described  by  the  following  set  of  equations: 


Here  e|j,  Sjj,  and  b|j  are  the  deviatoric  tensors  of  strain,  stress  and  backstress  respectively;  the 
difference  Sjj  -  b^j  is  often  referred  to  as  the  effective  stress.  The  quantities  I2,  J2>  and  K2  are  the 
second  order  invariants  of  the  deviatoric  tensors  of  strain,  dimensionless  effective  stress,  and 
strain  rate,  respectively;  these  quantities  are  defined  below: 

’  3  ®y  *  ^2  ~  2  ®y  ’  ^2  =  2  ®y 

“  ®y  *  3  ^y 


^  6jj*l  if  i=j,  6|j=0  if  i^j,  io=l  A3 
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Thus  the  growth  of  inelastic  strain  is  a  function  of  stress,  backstress,  and  strain  rate.  Also,  note 
here  that  plus  sign  appearing  with  the  radical  sign  of  the  square  root  of  the  invariants  in  Eqs.  (3)- 
(4)  indicates  that  the  square  root,  once  taken,  is  to  be  positive  (i.e.  the  absolute  value  of  the 
square  root).  The  tensor  e^j  is  known  as  the  distortional  component  of  strain  because,  by 
definition,  it  subtracts  the  dilatational  component  of  deformation  out  of  the  strain  tensor  Ejj . 
Therefore,  I2  represents  a  measure  of  volumetric  distortion  that  is  invariant  with  respect  to 
coordinate  transformations,  and  this  will  be  an  important  quantity  in  the  forthcoming  discussion. 
The  material  constants  in  Eqs.  (2)*(4)  are: 

E:  Young's  extensional  elastic  modulus 
v:  Poisson  ratio  of  elastic  material 

Of,:  minimum  axial  stress  necessary  to  activate  the  damping  mechanism 
a:  constant  which  determines  the  slope  of  the  inelastic  region 
=  Ey/(E  -  Ey),  where  Ey  is  the  inelastic  slope 

n:  constant  controlling  the  sharpness  of  transition  firom  elastic  to 
inelastic  behavior 

fp:  constant  controlling  the  size  of  the  hysteresis  loop 

a:  constant  controlling  the  amount  of  elastic  recovery  during  unloading 

Also,  Eqs.  (3)  and  (4)  contain  two  special  functions:  the  enor  function,  erf(  ),  and  the 
unit  step  function,  {u(  )}.  Simply  stated  the  purpose  of  the  error  and  unit  step  functions 
contained  in  Eq.  (4)  is  to  allow  for  the  recovery  of  accumulated  inelastic  strain  during  unloading, 
and  thus  simulate  the  unique  behavior  of  superelastic  materials  [13,14]. 

Let  us  take  a  moment  to  explain  the  role  of  the  inelastic  response  in  the  modeling  of 
strain  dependent  damping.  The  inelastic  component  of  strain  is  responsible  for  the  dissipation  of 
energy  that  takes  place  in  cyclic  loading.  Equations  (l)-(4)  have  been  used  to  represent  the 
macroscopic  stress-strain  behavior  of  shape  memory  alloys,  and  especially  superelastic  materials 
[13,14].  The  hysteretic  character  of  superelasticity  is  macroscopically  similar  to  that  of 
nonlinear  anelasticity  except  that  the  respective  stress  and  strain  levels  of  each  type  of  response 
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are  different  by  many  orders  of  magnitude.  Therefore,  the  inelastic  response  governed  by  Eqs. 
(3)-(4)  can  be  used  to  repi'esmt  the  maaoscopic  effect  of  a  nonlinear  anelastic  damping 
mechanism. 

By  using  Eqs.  (l)-(4),  a  number  of  special  cases  can  be  considered.  First  let  us  consider 
the  cases  of  uniaxial  tension'Compression  and  pure  shear  loading.  The  state  of  uniaxial  loading 
(superscript  u)  is  described  by: 
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Here  e,  a,  and  ^  are  the  axial  strain,  stress,  and  backstress  in  the  x  direction  of  Cartesian  space, 
respectively.  Also  the  lateral  strain  and  strain  rate  induced  by  the  Poisson  effect  (>pe  and  -pi:) 
are  associated  with  the  coefficients  p  and  p  respectively.  Strictly  speaking,  p  and  p  are  neither 
constant  nor  equal  due  to  the  uonlinear  effect  induced  by  the  damping  mechanism.  In  order  to 
evaluate  these  coefficients  the  lateral  strain  and  strain  rate  are  decomposed  into  elastic  and 
inelastic  parts  and  we  will  assume  that  volume  changes  induced  by  axial  loading  are  associated 
only  with  elastic  deformation.  Thus  the  inelastic  part  of  the  strain  and  strain  rate  are  associated 
with  incompressible  behavior  as  is  done  in  the  theories  of  plasticity  and  viscoplasticity  [17]. 

This  assumption  is  plausible  since  the  damping  mechanisms  involve  movement  of  dislocations, 
point  defects,  domain  walls,  or  polymer  chains  none  of  which  induce  a  change  in  volume.  Using 
this  assumption,  the  elastic  component  of  the  lateral  strain  is  related  to  the  elastic  axial  strain  by 
the  elastic  Poisson  ratio  v,  and  the  inelastic  component  is  related  to  the  axial  inelastic  strain  by 
the  Poisson  coefficient  of  incompressible  deformation  (which  is  0  J);  therefore  -pe  -  .5e“. 
Similarly,  the  lateral  strain  rate  is  -pe  — ve^^  •  .5e^.  Using  these  relations  it  can  be  shown 
(Appendix  A)  that  p  and  p  are: 
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Thus  (1  and  p  are  clearly  variable  coefilcients  which  are  not  necessarily  equal  to  each  other. 
For  the  conditions  of  pure  shear  loading  (superscript  s)  we  have: 
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Here  y  and  y  are  the  engineering  shear  strain  and  strain  rate,  r  is  the  shear  stress,  and  |  is  the 
shear  backstress  in  the  xy  plane  of  Cartesian  space. 

By  using  the  appropriate  stress,  backstress,  and  strain  tensors,  as  well  as  their  respective 
deviators  and  associated  invariants,  Eqs.  (l)-(4)  produce  the  following  uniaxial  equations: 


6  =  E 
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H  the  behavior  u  only  a  small  departure  from  elasticity  then  p.  “  p  »  v;  conversely  if  a  condition 
of  strain  and  strain  rate  exists  where  inelastic  behavior  dominates  and  where  o/e  «  E  and  do/de  « 
E  then  the  response  is  essentially  incompressible  with  p  »  p  *  .5.  Now,  since  p  and  p  appear 
only  in  terms  associated  with  inelastic  strain,  the  previous  equations  can  be  simplified  by  setting 
both  Poisson  coefficients  equal  to  0.5.  This  does  not  greatly  affect  the  numerical  results  since 
these  terms  will  be  significant  only  in  the  inelastic  region  (i.e.  when  oJE  <  e).  Thus  we  obtain: 


d  =  E 


£-1^1 


r  o 

P  =  Ea  j^E  -  g  +  fj  erf(aE)  {u(-ee)} 


(5) 

(6) 
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For  the  case  of  pure  shear  loading  the  governing  differential  equations  are: 


where 


t  =  G 


MyI 


{»>(‘YY)} 


G  =  rr - r  is  the  elastic  shear  modulus 

2(1  +v) 


(7) 

(8) 


is  the  shear  stress  whereupon  the  damping 
mechanism  is  activated 


Note  that  falls  out  of  the  formulation  automatically  in  a  manner  that  is  consistent  with  the 
theory  of  maximum  distortional  strain  energy  [18].  This  is  because  Eq.  (3)  is  dependent  on  the 
stress  gradient  of  a  potential  function  [14,16]  that  contains  a  Bingham  type  condition  for  the 
onset  of  the  inelastic  damping  mechanism. 

As  noted  before,  the  damping  mechinisms  become  saturated  at  an  axial  strain  of  Eg.  This 
differs  from  SMA  hysteretic  behavior  and  therefore  the  model  must  be  modified.  The 
modification  consists  of  specifying  that  the  growth  of  inelastic  strain  be  stopped  at  a  saturation 
limit  beyond  which  linear  elastic  behavior  once  again  takes  place.  To  do  this  the  second  term  in 
Eq.  (5)  is  multiplied  by  the  unit  step  function  {u(Eg  •  {eD}.  At  levels  where  e  <  Eg,  {u(Eg-|ED}=l 
and  inelastic  growth  of  strain  may  proceed,  but  when  e  a:  Eg,  {u(Eg-|ED}=:0  and  continued  loading 
beyond  Eg  is  elastic.  Similarly,  for  saturation  of  damping  mechanisms  in  shear  the  second  term 
on  the  right  hand  side  of  Eq.  (7)  is  multiplied  by  {u(Yo  -  tyD}.  where  Yq  is  the  shear  strain  of 
saturation. 

Since  we  are  modeling  materials  which  are  both  homogeneous  and  isotropic,  the 

saturation  process  in  more  general  three»dimensional  loadings  can  be  represented  by  a  value  of 
the  invariant  ^  which  is  determined  from  either  Eg  or  Yg'.  we  will  call  this  value  Thus 
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saturation  is  determined  by  a  specific  amount  of  distortion,  and  the  unit  step  functions  {u(e(,‘|eD} 
and  {u(Yo-|y|)}  can  be  obtained  from  {u(-^3ii^  -  yj^}. 

To  summarize,  the  equations  that  are  set  forth  to  represent  the  respective  three 
dimensional,  axial,  and  shear  responses  of  nonlinear  anelastic  materials  are  as  follows: 


3D: 

(5) 

bij  =  3  Eo  erf  (ja  (u  (-I2) ) 

(10) 

Axial: 

o  =  E  £-|e|  {u(Eo-N)} 

c  ^  c  ^ 

(11) 

P  =  Ea  E-^+fj erf(aE) {u(-ee)} 

(12) 

Shear: 

X  -  H  fT  -  £1 

t  =  G  y-IyI  ,  {u(yo-IyI)} 

(13) 

^  =  T  ^  {u(-YY)} 

(14) 

where 

2(1 +v) 

Yo  -  **  Strain  limiting  inelastic  growth  in  shear 

The  hysteresis  loops  associated  with  axial  and  shear  behavior  are  produced  by  numerical 
integration  of  Eqs.  (11)-(12)  and  (13)-(14),  and  by  specifying  a  sinusoidal  history  of  strain  input 
with  amplitudes  of  Ep  in  the  axial  case  and  Yp  in  the  shear  case  (Appendices  B  and  C).  These 

amplitudes  were  specified  to  be  greater  than  the  saturation  strains  so  that  the  full  character  of  the 
predicted  response  could  be  plotted.  The  results  of  calculations  for  the  axial  and  shear  loading 
conditions  are  given  in  Figs.  4  and  5  respectively.  Both  figures  possess  the  same  characteristics: 
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elastic  behavior  dominates  in  the  region  of  the  origin  as  well  as  outside  the  limiting  strain,  and  a 

hysteresis  loop  is  manifested  in  the  full  cycle  of  strain  application.  The  area  enclosed  by  the 

hysteresis  loop  represents  the  energy  absorbed  by  the  material  per  cycle  of  oscillation,  and  this 

> 

quantity  is  denoted  as  AW.  An  elastic  modulus  of  E  =  28.5x10^  psi  and  saturation  strain  of  Eg  = 
0.0001  were  selected  based  on  the  elastic  modulus  and  approximate  strain  of  peak  damping  in 
Fe-Cr  alloys  (see  Hg.  1).  The  remaining  material  constants  used  in  the  calculations  which 
generp.ted  Figs.  4  and  5  were  not  selected  to  reproduce  the  behavior  of  any  specified  damping 
material;  rather  they  were  selected  to  approximate  a  typical  pattern  of  hysteresis  in  nonlinear 
damping  materials  and  to  permit  investigative  analyses.  It  should  be  noted  that  the  numerical 
results  generated  by  the  constitutive  equations  were  obtained  by  Runge-Kutta  fourth  order 
integration,  and  these  results  for  the  material  response  were  strain  rate  independent.  Also  the 
elastic  and  inelastic  material  properties  (i.e.  E,  Ey,  and  Oj)  are  accurately  reproduced  in 

numerical  calculations  [14]. 

By  having  numerical  results  of  the  type  just  presented,  we  can  now  compute  the  damping 
according  to  the  definition  of  the  loss  factor,  t],  defined  as 

AW 
^  ~  2jiW 


where  AW  is  defined  as  given  above,  and  W  is  a  measure  of  stored  energy  most  often  selected  as 


max 


(16) 


By  using  Eqs.  (11)-(12)  and  (13)-(14)  in  calculations  for  the  cyclic  material  response  over  a 
range  of  peak  axial  and  shear  strains,  and  computing  the  loss  factor  associated  with  each  peak 
strain  according  to  Eqs.  (15)  and  (16),  the  general  character  of  the  damping  vs.  strain  diagram  of 
nonlinear  materials  was  produced  (computer  algorithms  are  given  in  Appendices  D  and  E);  this 
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is  shown  in  Fig.  6.  Note  that  both  curves  possess  the  characteristic  damping  peak  associated 
with  nonlinear  damping  materials. 

However  the  separate  curves  in  i^g.  6  representing  axial  and  shear  loading  differ 
significantly  with  respect  to  one  another,  and  this  can  make  the  results  difficult  to  interpret.  In 
an  efibrt  to  understand  these  results,  the  amount  of  energy  absorbed  in  each  loading 
configuration  was  evaluated.  The  amount  of  energy  absorbed  is  plotted  as  a  function  of  peak 
strain  in  Fig.  7.  Above  the  respective  axial  and  shear  saturation  strains  the  amount  of  energy 
absorbed  by  the  material  is  essentially  the  same  for  both  loading  configurations,  their  difference 
being  less  than  1%.  Therefore  the  difference  in  the  character  of  the  two  separate  responses  must 
be  due  to  other  factors. 

It  turns  out  that  the  plots  given  in  both  Figs.  6  and  7  are  misleading  because  the  abscissa 
of  these  figures  represents  values  of  strain  associated  with  separate  axial  and  shear  loading 
conditions,  and  the  strains  associated  with  these  separate  conditions  are  not  equivalent  Since 
damping  is  a  material  property,  it  is  desirable  to  plot  strain  dependent  material  damping  values 
using  a  measure  that  will  unify  the  curves  from  separate  tests.  Such  a  method  would  also 
provide  a  basis  by  which  to  present  and  compare  damping  data  for  a  variety  of  materials. 
Therefore  another  measure  of  deformation  equivalent  to  both  types  of  loading  needs  to  be 
employed.  One  such  possibii.'ty  is  to  use  a  measure  of  distortion  rather  than  strain.  Let  us  define 
an  equivalent  strain  I  as  follows: 

6=-v/35  (17) 

This  measure  is  similar  to  the  effective  plastic  strain  in  plastically  deforming  materials  [17].  It  is 
clear  that  e  has  a  physical  meaning  that  is  independent  of  the  choice  of  coordinate  axes  since  it  is 
based  on  the  invariant  fy,  which  is  the  second  invariant  of  the  deviatoric  strain  ey.  Therefore  e  is 

an  invariant  measure  of  distortion. 
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By  considering  the  separate  conditions  of  axial  and  shear  loading,  and  by  taking  the 
variable  Poisson  coefficient  to  be  the  elastic  constant  v  for  axial  case,  the  equivalent  strains  for 
each  condition  are: 

E”  s  (1  +  v)  E  (uniaxial  loading) 

6®  =  (shear  loading) 

Using  the  peak  equivalent  strains  of  axial  and  shear  loading  in  place  of  the  peak  strains  used  in 
Figs.  6  and  7  produces  a  more  consistent  pattern  of  results.  This  is  shown  in  Fig.  8  where  the 
energy  absorbed  as  a  function  of  £  is  in  very  good  agreement  along  the  entire  abscissa  for  both 
cases.  In  Fig.  9  the  loss  factors  of  the  axial  and  shear  loading  cases  are  also  plotted  against  the 
peak  equivalent  strain.  Even  though  the  peaks  of  the  separate  curves  are  not  of  equal  magnitude 
the  results  in  this  figure  are  now  very  similar;  indeed  the  rise  and  fell  of  each  curve  follow  the 
same  trend  and  the  peak  of  each  damping  curve  occurs  at  approximately  the  same  level  of 
distortion. 

The  peak  of  the  loss  factor  curve  for  shear  loading  is  slightly  higher  than  the  peak  of  the 
axial  loss  factor  curve  due  to  a  smaller  value  of  stored  energy  in  the  shear  loading  case.  This  is 
another  complication  that  arises  from  the  nonlinear  nature  of  the  of  the  stress-strain  material 
response  and  (for  equal  amounts  of  distortion)  causes  the  value  of  peak  shear  stress  to  be  lower 
than  the  peak  axial  stress.  Consequently,  at  equal  levels  of  distortion,  the  measure  of  stored 
energy  (Eq.  (17))  will  be  larger  in  axial  loading  than  in  shear  and  this  will  cause  the  loss 
modulus  in  shear  to  be  greater  than  the  1(^  modulus  in  axial  loading. 

Bending  and  torsion  are  common  configurations  in  which  to  measure  damping.  The 
solid  beam  has  length  L  and  rectangular  cross-section  of  width  b  and  thickness  h.  Although  the 
stress-strain  re^nse  is  nonlinear  we  can  consider  both  cases  in  a  simple  fashion  because  the 
response  takes  place  in  a  manner  which  gives  symmetric  behavior  for  positive  and  negative 
strains.  When  considering  bending  and  torsion  problems  with  pronounced  plastic  deformation 
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and  noDsymmetrical  stress-strain  responses,  then  special  considerations  must  be  made  when 
computing  the  acting  moments  [17], 

Schematic  illustrations  of  bending  and  torsion  are  shown  in  Fig.  10.  The  torsional  solid 
shaft  has  length  L  and  circular  cross-section  of  radius  R.  Note  that  the  strain  profiles  in  each 
geometry  are  linear,  passing  through  zero  at  the  position  of  the  neutral  axis  of  the  beam  and 
starting  at  zero  at  the  center  of  the  shaft  Also  note  that  Ep  is  the  value  of  the  axial  strain  at  the 
beam  surface  while  Yp  is  the  value  of  the  engineering  shear  strain  at  the  shaft  surface.  The 
angles  6  and  ^  are  the  curvature  of  the  bending  beam  and  the  angle  of  twist  of  the  shaft 
respectively.  Because  the  problems  under  consideration  involve  only  small  strain,  the  following 
simple  relations  can  be  used  to  compute  the  moments  and  angular  displacements  for  the  beam 
and  shaft  geometries  respectively: 


M  =  -J  y  odA 
A 


and 


(beam) 


T  =  J  rt  dA  and 

A 


(shaft) 


Here  y  is  the  vertical  distance  from  the  neutral  axis  of  the  beam  cross-section,  o  is  the  axial 
stress  in  the  longitudinal  fibers  of  the  beam,  and  M  is  the  resultant  moment  bending  the  beam; 
for  the  shaft  r  is  the  distance  from  the  center  of  the  circular  cross-section,  x  is  the  shear  stress 
due  to  torsion,  and  T  is  the  resultant  torque  twisting  the  shaft.  The  equations  given  above 
relating  6  and  (j>  to  sample  dimensions  and  sur&ce  strain  amplitude  are  easily  deducible  from 
simple  geometrical  arguments  which  involve  the  knowledge  of  a  linear  strain  profile  and  the 
restriction  of  small  strains. 

In  the  numerical  analyses  associated  with  these  bending  and  torsion  tests,  the  surface 
strains  of  each  test  sample  were  specified  to  act  sinusoidally  in  time.  Each  geometry  was 
subdivided  into  a  large  number  of  finite,  but  thin,  subsections  (i.e.  the  infinitesimal  distances  dy 
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and  dr  in  Fig.  10  were  replaced  by  small  but  finite  distance  Ay  and  Ar  respectively).  Also,  the 
strain  distribution  for  each  finite  subsection  was  assumed  to  be  constant  over  the  subsection 
thickness  and  the  value  of  the  strain  was  taken  as  the  value  of  the  strain  profile  at  the  center  of 
the  subsection.  Knowing  the  strain  profile  of  the  cross-section  of  each  geometry  and  the  history 
of  the  respective  surface  strains,  the  stress  history  for  each  subsection  of  the  geometry  was 
computed  numerically.  Specifically,  Eqs.  (11)-(12)  were  integrated  to  give  the  stress  profile 
time  history  of  the  bending  beam  and  (13)-(14)  were  integrated  for  the  shear  stress  profile  time 
history  of  the  shaft  (see  Appendices  F  and  G).  Then  the  following  formulas  were  used  to 
compute  the  resultant  moment  and  torque  histories  of  the  beam  and  shaft; 

N 

M  =  -b  2)  Yi  Oi  Ay 

i=l 

N 

T  =  2ji2)(ri)2xiAr 


where  N  is  the  number  of  subdivisions  making  up  the  cross-sectional  geometry  and  where  the 
subscript  i  indicates  reference  to  the  location  of  a  single  subsection. 

The  loss  factor  of  each  sample  geometry  was  then  calculated  for  a  specified  value  of 
surface  strain  amplitude  according  to  Eq.  (15)  where  AW  was  determined  by  the  area  enclosed 
by  the  resultant  moment  vs.  angular  displacement  hysteretic  response  and  W  was  determined  by 


(beam) 


(shaft) 


The  damping  values  which  were  computed  in  this  way  were  found  to  be  independent  of  sample 
geometry,  i.e.  for  a  given  surface  amplitude  the  ratio  of  AW  to  W  remained  constant  for  changes 
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in  cross-sectional  size,  sample  length  or  both.  This  effect  is  due  to  the  fact  that  a  proportional, 
volume  of  material  is  undergoing  deformation  wherein  the  damping  mechanism  is  activated  and 
this  proportion  of  volume  is  constant  irrespective  of  sample  size  for  both  simple  bending  and 
torsion. 

In  Fig.  11  the  loss  factor  was  plotted  against  the  sinface  amplitude  for  both  the  bending 
and  torsion  cases  by  repeating  the  calculations  over  a  range  of  surface  amplitudes.  Note  that  the 
character  of  the  damping  vs.  surface  amplitude  curves  are  different  with  respect  to  one  another. 
This  is  analogous  the  trend  shown  earlier  in  Fig.  6  for  one  dimensional  behavior.  Also,  by 
comparing  Fig.  11  to  Fig.  6  it  is  clear  that  the  character  of  the  damping  vs.  peak  strain  curve  of 
each  sample  is  quite  different  than  that  corresponding  to  the  respective  one-dimensional  material 
point  responses,  lliis  is  due  to  the  strain  dependent  nature  of  the  damping  and  the  fact  that  strain 
is  nonuniformly  distributed  throughout  the  sample;  therefore  some  regions  of  the  geometry  may 
be  contributing  significantly  to  the  overall  damping  of  the  solid  sample  while  others  are  not. 

Using  Eq.  (11)  to  calculate  the  amounts  of  peak  equivalent  strain  at  the  surface  of  the 
bending  and  torsion  samples,  Fig.  12  shows  that  the  use  of  peak  equivalent  strain  gives  an 
improved  measure  of  correlation  in  the  same  manner  that  was  exhibited  earlier  for  the  one¬ 
dimensional  cases. 

Thus  presentation  of  nonlinear  damping  data  as  a  function  of  equivalent  strain  rather  than 
as  a  function  of  sample  strain  can  be  very  useful.  It  is  probably  most  useful  in  comparing 
damping  data  obtained  by  different  test  methods.  It  may  also  be  useful  in  design  work  where  the 
dynamic  strains  in  a  vibrating  part  or  member  are  known.  Along  this  line,  let  us  briefly  consider 
an  example  where  a  designer  wishes  to  use  a  high  damping,  but  nonlinear,  material  in  an 
application  where  bending  is  the  primary  mode  of  deformation,  and  suppose  that  damping  data  is 
available  only  from  torsional  tests.  If  the  vibrational  strain  levels  to  be  expected  in  service  can 
be  deduced  from  load  and  design  analyses,  then  these  strain  levels  can  be  converted  to  the 
measure  of  equivalent  strain  introduced  in  this  paper.  The  designer  would  then  be  able  to 
estimate  whether  or  not  the  material  damping  will  be  in  a  range  of  peak  performance  for  the 
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application  of  interest  by  applying  the  same  conversion  to  the  peak  shear  strains  of  the  torsional 
damping  data. 

SUMMARY 

The  work  presented  in  this  paper  includes  three  major  aspects;  1)  modeling  of  nonlinear 
(or  strain  dependent)  damping  behavior  via  constitutive  equations,  2)  relating  the  damping  of  the 
material  to  the  damping  of  a  test  specimen,  and  3)  a  way  of  improving  correlation  of  nonlinear 
damping  data  via  use  of  equivalent  measures  of  distortion.  These  efforts  were  conducted  in 
order  to  gain  a  better  understanding  of  macroscopic  nonlinear  high  damping  material  behavior 
and  also  to  obtain  a  means  in  which  to  better  correlate  damping  data  from  tests  which  use 
different  sample  geometries.  The  modeling  scheme  applies  to  homogeneous  isotropic  materials 
and  is  adapted  from  a  constitutive  model  of  the  viscoplastic  type  through  incorporation  of 
constants  that  represent  the  onset  of  damping  mechanisms.  Also  the  model  was  modified  to 
include  damping  mechanisms  that  become  saturated  after  a  given  amount  of  strain.  Analyses 
were  made  to  calculate  the  loss  factor  of  the  common  damping  test  configurations  of  bending 
and  torsion.  To  do  this  material  point  relationships  were  used  at  a  large  number  of  points 
making  up  the  cross-sectional  geometry.  In  this  way  it  was  possible  to  relate  the  damping  of  the 
material  to  the  damping  of  the  specimen.  Hie  results  did  not  depend  on  the  relative  dimensions 
of  the  sample  geometry;  rather  the  calculated  loss  factors  depended  only  on  the  mode  of 
deformation.  The  results  showed  that  the  strain  dependent  damping  associated  with  each  test 
were  difficult  to  compare  when  plotted  solely  against  the  peak  surface  strain  of  the  sample 
geometry.  This  is  because  the  peak  strains  that  correspond  to  each  of  these  test  configurations, 
namely  axial  and  shear  strain,  are  different  from  one  another.  However  if  an  invariant  measure 
of  peak  sample  distortion  is  used  in  place  of  peak  sample  strain,  then  the  correlation  of  the 
nonlinear  damping  of  separate  bending  and  torsion  samples  improves  considerably.  Such  an 
improved  capacity  for  the  correlation  of  nonlinear  damping  data  is  very  useful  for  comparison  of 
data  obtained  from  different  tests.  Future  research  will  include  the  modeling  of  specific 
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nonlinear  damping  data.  Also,  constitutive  model  material  parameters  that  are  physically 
motivated  by  the  microstructure  will  be  studied. 
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LOSS  FACTOR 


Fig.  1.  Strain  Amplitude  Dependent  Damping  in  Fe-Cr  Based  High  Damping  Alloys. 


Cyclic  ff-e  with  a  Saturated  Damping  Mechanism  Outside  £o 


Strain 


Fig.  2.  Generalized  Macroscopic  Hysteresis  of  Nonlinear  Damping  Materials. 
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LOSS  FSCTOS  VS 


Fig  3.  Strain  Dependent  Damping  of  Qi-Mn  in  Separate  Bending  and  Torsion  Tests  [3]. 


Hysteretic  Model  of  Axial  Loading  with  Elasticity  Past  £o 
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Axial  Strain  (in/in) 

Fig.  4.  Hysteretic  Behavior  Calculated  for  Pure  Axial  Loading. 
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Shear  Stress  (psi) 
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Energy  Absorbed  by  Axial  and  Shear  Hysteretic  Models 


Peak  Strain  (Axial,  Shear) 

Fig.  7.  Energy  Absorbed  in  Axial  and  Shear  vs.  Peak  Axial  and  Shear  Strains. 
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a=80000  fT=4xi0"* 

i.l  J.i.1  I  :  M  M  I  I  1  t  H  I  I  I  J_1  I  1  .1  1  I  1  1  >!  1  I  I  1  I  I  1  I  I  I  I  I  1  I  I  I  I  I  1  I . Ill' 


0.0000  0.0001  0,0002  0.0003 

Peak  Equivalent  Strain  (Axial,  Shear) 

Fig.  8.  Energy  Absorbed  in  Axial  and  Shear  vs.  Peak  Equivalent  Strain. 
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actor 


Peak  Equivalent  Strain 


Fig.  9.  Amplitude  Dependent  Damping  in  Axial  and  Shear  vs.  Peak  Equivalent  Strain. 
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Bending  Beam  (Length:  L,  Width:  b.  Thickness:  h) 


M  =  -  J  y  o  dA  -  -b  2  yj  Oj  Ay 


Twisting  Shaft  (Length:  L,  Radius  R) 


A 


T=JrxdA-2iiy  (rj)2Ti  Ar 


Fig.  ly.  Schematic  Dravnng  of  Strain  Profile  in  Bending  and  Torsion  Geometries. 
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strain  Dependent  Damping  in  Beam  and  Shaft  Samples 


Peak  Strain  at  Specimen  Surface  (in/in) 

Fig.  11.  Amplitude  Dependent  Damping  in  Bending  and  Torsion  vs.  Peak  Surface  Strain. 


Strain  Dependent  Damping  in  Beam  and  Shaft  Samples 


Peak  Equivalent  Strain 


Fig.  12.  Amplitude  Dependent  Damping  in  Bending  and  Torsion  vs.  Peak  Equivalent  Strain. 
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APPENDIX  A 


In  this  appendix,  we  will  examine  the  Poisson  effect  associated  with  isotropic  material 
deformation  with  nonlinear  inelastic  effects.  To  do  this  first  consider  the  strain  tensor  e^j 

associated  with  pure  axial  stressing,  i.e.  without  off-diagonal  (shear)  strains.  The  total  strain  is 
made  up  of  separate  elastic  and  inelastic  parts  as  follows: 


where  the  elastic  part  (superscript  el)  is  related  to  the  strecs  by  elasticity  theory  and  where  the 
inelastic  part  is  associated  with  incompressible  deformation  mechanisms  (i.e.  mechanisms  which 
operate  without  any  associated  volume  changes;  e.g.  dislocation  glide).  Therefore,  for  the  case 
of  uniaxial  (superscript  u)  stressing 

re«l  0  0  1  6“  0  0 

Eij  =  =  0  -vecI  0  +  0  -|e“  0 

ij  ij  ^ 

.  0  0  -VE«1J  0  0 

Note  that  the  total  axial  strain  is  specified  to  be  aligned  along  the  x  coordinate  of  Cartesian  space 
with  E®1  =  o/E.  This  in  turn  leads  to 

ecI  4.  gin  0  0 

0  -ve®^-|e“  0 

0  0  -VE®^  - 

Thus  it  is  clear  that  the  Poisson  induced  strains  are  composed  of  an  elastic  part  which  is 
related  to  the  elastic  strain  by  Poisson's  ratio  and  an  inelastic  part  which  is  related  to  the  inelastic 
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strain  by  the  coefficient  1/2.  By  this  process  it  is  clear  that  the  inelastic  part  of  strain  satisfies 
the  condition  of  incompressibility. 

Now,  by  using  =  o/E  and  e  =  +  £“  let  us  evaluate  the  Poisson  induced  strains  in 

more  detail.  First  note  that 


Thus,  by  defining  a  variable  Poisson  coefficient  p.  as 


1  1  fl  1  o 

^  =  2-e[2-^U 


it  follows  that 


-pe  =  -  ve®l  - 


(Al) 


(A2) 


u 

®ij 


E  0  0  • 

0  -pE  0 
.0  0  -pE. 


(A3) 


.u 


Now,  from  Eq  (A3)  we  can  compute  the  strain  rate  tensor  e  ■ 


E 


u 

ij 


E  0  0 

0  -ps  -  pE  0 

0  0  -pE  -  pE  , 


From  this  let  us  proceed  to  evaluate  the  Poisson  induced  strain  rates  in  more  detail.  Considering 
Eq.  (Al)  it  follows  that 
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1  r  1  1  r  o  .  1.1 


and  this  leads  to 


-|1£  -  |XE 


=  -26  + 


=  -veel-|(  E-eel)  =  .vEel-^Ein 


(A4) 


Thus  the  Poisson  induced  strain  rate  decomposes  into  two  parts:  one  part  is  an  elastic  component 
and  is  related  to  the  axial  strain  rate  by  the  Poisson  ratio,  and  an  inelastic  part  that  is  associated 
with  the  inelastic  strain  rate  in  a  manner  consistent  with  incompressibility. 

By  using  e®!  =  6/E  and  e  =  e®!  +  e“  Eq.  (A4)  becomes 


Thus  the  variable  Poisson  coefficient  associated  with  strain  rate  (which  is  denoted  here  as  p)  is 
as  follows: 


1  1  rl  1  do 

P"2’e[2‘''J  dE 
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and 


E 


u 

ij 


E  0  0  ' 

0  -PE  0 
0  0  -PE. 


for  axial  stressing. 
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APPENDIX  B 


FORTRAN  MAIN  PROGRAM  (all  subroutines  and  function  subprograms  that  are  not 
specifically  related  to  the  calculations  being  made  here  are  given  in  Appendix  H): 


■»»">  Program  Hanet  UNIAXIAL. FOR 

■>  This  program  will  carry  out  a  Runge-Kutta  integratiou  on  the 
stress-strain  equations  of  the  SMA  hysteretic  model  for  eocial 
loading.  Note  that  inelastic  behavior  is  suppressed  beyond  EFSLNO. 
The  strain  history  is  prescribed  (i.e.  input)  as  sinusoidal. 

Stress  is  the  output  variable.  Strain  and  stress  pairs  (i.e.  the 

hysteresis  curve)  are  written  to  data  files  for  plotting 

with  GRAPHER  (which  reads  ASCII  data  arranged  in  column  pairs). 


NOMENCLATURE  t 


REAL  CONSTANTS! 

A  Material  constant  controlling  shape  of  hysteresis 

ALFA  “Ey/(E-Ey)  where  Ey  is  the  inelastic  modulus 

AMP  Amplitude  of  strain  input 

E  Young's  modulus 

EPSLNO  Strain  beyond  which  inelasric  growth  is  supro&sed 

FREQ  Frequency  of  cyclic  strain  input 

FT  Material  constant  controlling  size  of  hysteresis 

N  Overstress  power  (controls  sharpness  of  transition  to  inel.) 

VNEW  Poisson's  ratio 

VNENIH  Inelatlc  Poisson  coefficient 

Y  Stress  where  damping  mechanisms  are  activated 

INTEGER  CONSTANTS! 

NCiCLE  No.  of  cycles  of  oscillation 

:<PPCYC  No.  of  points  per  cycle  to  be  used  in  integration 
CHARACTER  STRINGS! 

FILENAME  String  for  filename  asnigruaent  to  a  FORTRAN  unit  number 

SALF  string  for  material  constant  ALFA 

SA  string  for  material  constant  A 

SE  string  for  Young's  modulus  E 

SEO  string  for  limiting  strain  EPSLNO 

SEP  string  for  pea)c  strain  AMP 

SFT  string  for  material  constant  FT 

SN  string  for  overstress  power  N 

SY  string  for  stress  where  damping  mechanisms  activate:  Y 

TITLE  descriptive  title  for  ru.^ 

VARIABLES! 

ENGABS  Energy  2d}sorbed  per  cycle  of  oscillation 

ENGSTO  Energy  stored  (>.5  *  max.  strain  ^  max.  stress) 

ETA  Loss  factor 

SDC  Specific  damping  capacity 


VARIABLE  ARRAYS! 

EPSLOH(K)  Axis  Strain  (index  K  rep.  time) 

STRESS (K)  Axia..  stress  (index  X  rap.  time) 

Z(X)  Axial  stress  passed  frocn  INTFUN  (index  R  rep. 


time) 


SUBROUTINES! 

DERIV 

INTFUN 

MAXIM 

STRLEN 


Contains  the  differential  eqs.  (invoked  by  INTFUN) 
Integration  routine  (4th  order  forward  Runge-Kutta) 
Determines  the  min.  t  suuc.  values  of  an  array 
Counts  the  numbar  of  characters  in  a  string 


FUNCTION  SUBPROGRAMS! 

FRP(X)  Error  function  of  X 

FACT(K)  K  factorial 

SGN(V)  Signum  function  of  V 

UNIT(X)  Unit  step  function  of  X 


IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  NEQ 

REAL*8  DERIV,FLOAT,T,TEND,TOL,Z(20),H 
DIMENSION  EPSLON(1001),STRESS(1001) 
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CHARACTZR*20 

CHaRSiCT3R*70  TS7L2 ,  SY , SB , SH, SALP , SEO , fi&, SPT,  32P 
CH>ilACn!R*(*)  SARi,SAR2,SAR3,aaR4,SRR5,SAR6,SAE7,SASP 
PARAMETER  (SARl-'eO~’) 

PARAMETER  {S&R2-'Y»’) 

PARAMETEC  (SAR3-*E-') 

PARAMETER  { SAR4« ' n- ' ) 

PARAMETER  (SAR5«’aifa»' ) 

PARAMETER  { SAR6» ' a- ' ) 

PARAMETER  (SAa7-’fT-’) 

PAHAMEfER  ( SARP- *  tip- * } 

COKMOH/BiiOKl/AHP ,  OffiGA ,  Y ,  E ,  VHEWIB  ,  AUA,  H ,  E?sr  J!0 ,  A,  PT 


HEQ  »  1 
C 

T  ‘T  <j.O 
I(i>  -  0.0 
C 

Intaractj.va  input  of  fllencnes  for  the  material  data  input  file 
C  iu2d  for  output  files. 

C 

WRITE (»,*', 'Enter  the  name  of  your  input  data  file.' 

IPL«9 

1  Wai)(»,2,P3XEMAME 

2  PORMAT(A) 

WRITE  ( ‘ , 3 ) PILEHAKE 
3FORMAT{‘  ',3X,A20) 

I?(IFL.EQ.9)THEd 

OPEN  ( IPL ,  I’jXE-FXIiEBAKE ,  STATUS-  *  CHKaOHK ' ) 

ELSE 

OPEN ( IPL, PILE-PILENAME , STATUS- ‘ OHEHOWN ' ) 

ENDIP 

IP  (IPL. EC. 9)  THEN 

WSi'rE(»,») 'Enter  acae  for  the  plot  file  for  STRESS  vb.  EPSLON' 
IPL-iO 
O  )  TO  1 
Ewnip 

If  (I?L.E0,10)  then 

WRITE C*,*)  ‘Enter  naae  for  the  aumiary  text  file’ 

1PI/-11 
GO  TO  1 
ENDIP 
REWIND  09 
REWIND  10 
REWIND  11 
C 

C— »>>  Read  input  quantities  froai  input  file. 

C 


READ  (9. 5)  TITLE 
REAi:(0,‘{A)')SY 
READ(9,  '(A)  '  )S’-i 
RRAD(9, ■ (A)  '  )SH 
RE/D{9, ' (A) ' jSALF 
S  FOR'4AT(A) 
read(9,*;e,vnew 
READ(9,*}Y 
READ(9,*)ALPA 
SEAD(9, »)S 
C 

>>  Interactive  input  of  other  naterial  parajseters,  strain  aaplituda, 
C  ano  nusber  of  cyciec  of  loading  to  be  used  in  calculations. 


WR1TE(*,*)'  Enter  EPSLHO' 

P£AD(5,^)EPS2iN0 

WRITE (•,»)'  Enter  character  string  for  EPSLNO  * 
READfS, '(A)  =  fSEO 
WRITE'*,*)  Enter  A  AND  FT* 

READ(5,*)A,FT 

WRITE (*,*)'  Enter  cDeracter  string  for  a' 
SSAB(5,  ■  (A)  •  ;-SA 


WRITE (*.*)■  Enter  character  string  for  f T ' 

REA0(5, • (A) ' )SPT 

WRITEC*,®;'  Enter  FREQ  and  AMP' 

RKAD{5,*)PREQ,.\MP 

WRITE {«,*)’  Ev.ter  character  string  for  Peak  Strain  (i.e.  AMP)' 
READ(5, • (A) ' )SRP 

WHITE (»,*)'  Enter  Kitaaber  of  ryclss  and  Susber  of  Point*  per  Cyclo' 
BEAD ( 5 , * ) SCYCLE , HPPC VC 
C 

C«”— — >  uetemljse  lecgth  of  various  character  strings  for  later  use. 

C 
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CAIX  ST:iIiEK(SBO,IBSBO,ZESEO) 

CALL  S7RLEN(SEP,IBSEP,IE5BP) 

CALL  STRLBH(SY,IBSY,IESY) 

CALL  STRLEN(SE,IBSE,XF.SB) 

CALL  STRLEH(SN,IBSN,1ESK1 
CALL  STKI2Ijf(SALP,IBSALP,IBSAL?) 

CALL  STRLEN(SA,ZBSA,1BSA) 

CALL  STRLEH(SPT,1BSFT,IESPT) 

CALL  S;5tLEH(SAP.l,ZBSARl,XS&AP.l) 

CALL  STItL2H(SAR2,IDSAR2,tBSAR2) 

CALL  S7BLBK(SAR3,IBSAR3,I£SAP.3) 

CALL  STRLEH(SAR4,IBSAR4,XESAR4) 

CALL  3T»LEK(SAR5,IBSAR5,IBSAR5) 

CALL  S'rBLBH(SAR«,XBSAR4,ZSSAR6) 

CALL  S7RLEN(SAR7,ZBSAB7,I£SAR7) 

CALL  STia,EH(SARP,ZBSAKP,IESARP} 

C 

ISDKl  -  rjSAAl  +  L2SE0  +  2  +  lESAPJ  +  lESEP 

IS0K2  “  IESAR2  +  lESY  +  2  +  IESAR3  +  lESB 

IS0M3  «  I£«iAR4  +  lECH  +  2  +  IESAR5  +  lESALP 

IS0M4  -  IFSAR6  f  IZhA  +  2  +  IESAR7  +  IBSPT 

C 

C>n»>»»>  Calculate  the  quantity  Pi«3. 14159267...,  (uid  other  para®etar8 


PI«DACOS(-1.9DOO) 

C 

OMBGA  -  2.»PT»FRBQ 

VH2WIH  -  .5  -  ALFA*(  .5-VKEW)/(l,+ALPA) 

SSTEPS  »  HCVCLE^KPPCre 
PRIHT  HSTEPS-' ,KSTEPS 
PERIOD  -  l./FREQ 
DEM  -  PERIOD/HPPCYC 
PRIHT  DELY“‘,DELT 
C 

C— >««.>  Set  Initial  coaditlonL. 
c 

STRESS(l)  -  0. 

EPSL0K{1)  >  0. 

5«aTE(*,S599S) 

99999  P0R«AT(4X, 'ISTEP',SX, •STRAIK',9X, *SrRELS’) 

C 

€•»•«'•»«>  Carry  out  the  nuB)-'3riL<il  integration 

c 

DO  10  ISTEP  -  l.KSTBPS 

CALL  mTFU<l(Z,T,rELT,HEQ) 

STKESd(ISTS?+i)  -  Z(l) 

BPSLOH(TdTEP+X)  ->  AMF  <»  DSIH<OKEGA«T) 

WRITS(»,  ■  (3X,I6,4(3X,S10.4)  )')TSTEP,EPSLOHJISTEP+l),r.a) 
10  CONTIKUP 


C 

C- 

C 

C 


“““>  Ccspute  dasping. 

CALL  MAXIH{ STRESS, SSTEPS, STRMIS.STEMAX) 

ENGABS  >  0. 

XSTART  »  KFPCYC/4.  +  1 
lEND  •»  XSTART  +  NPl>CYC-i 
DO  SO  I-ISTAKT, lEHD 

EHCABS  •>  ESGABS  .5  •  (STRESS(1+1)  +  STP.ESS(I))  * 

>  (EPSLOS{X+i5  -  ErSLOS(I)) 

50  COSTIHira 

EHGSTO  »  .5*STRMAX*AHP 
ETA  «•  ESGABS/ ( 2. <PI*EHGSTO) 

SDC  -  ETA* 2 •PI 

“>  Krite  resuito  to  output  data  file  (unit  10)  and  to  output 
tn'jTL  file  (unit  11). 


vniTEilO, •{lS,2{E16.e, ' 


’)}')( EPSLOH ( I ) . STRESS ( T ) , J”1 , hSTEPS ) 


fRlTEfil,  (5X, ‘*0  “,12, . ,A,A,2X,A,A  )*) 

>  ISUMl , SARI , SEO ( IBSBO : lESEO ) , SART, SEP/ XBSEP  J lESEP ) 

KRITE(n,‘(5X.  "I  •■,12, . ,A,A,2X,A,A,  ••  •'•)') 

>  ISUM2 , SAR2 , SY ( IBSY : XESY } , ^ AR3 , SE ( TBSE  J lESE ) 

HRlTE.ll, -(SE, • *2  ’■,I2,''  - • ’ , A , A, 2X,A,A, ’ ‘ * J ' ) 

>  XSaH3  ,  SFB4 ,  SN  ( IB5H I  IEs,S  >  ,  SAR5 ,  SALF  ( IBSALF I I2SALF ) 

WRITE ( 11 .  ■  J  f,S,  '  ■  3  ■  * ,  i: , . , A,A,  2X,A,  A,  ••-■•)• ) 

>  I3UH4  ,  SARC ,  SA(  IBSA:  lESA)  ,  SAR7  ,  SFT(  X8SFT I  lESPT) 

WRITS(11, ‘{SX,  ' ‘4  2S  ■•h-‘',L9.4,2X,  "SDC-’ ',29.4,  •*■'•')  •)ETA,SDC 
PRISS*,'  XjOSS  FACT0K-’ -ETA,  •  SDC-', SDC 


DTRC^ME-91/34 


lOOC  FORHa?(t70) 
lOOi  PORMAT(2A60) 

1200FORMAr(’  ‘.DXf 'I’.BS,!?) 

■,3X, •2',5X,I7) 
■,5X,Ei0.4,5S,E10.«) 
•,5X,K10.4,5S,Ei0.4,SX,E10.4) 


1201  miMAT(  ' 

1500  FOaMAT(  * 

1501  pomwri ' 
STOP 
END 


SOSROUIIKE  BERIV(£,T,ZtX)T> 

IH^LICIT  REAi*8  (A-'J,0-Z> 

X^TEGER  N2Q 

PJSA1*8  05ffiaA,‘r-l(20),Z50T(205,aH,XOJ«)0 
COHKON/BLOKl/AkP,OMEGA, Y, E . VEEWXK, ALFA.BH, EPSLHO, A, PT 
C 

EPSIiOa  “  AKP  »  CSIS(0}IEGA»T) 

EPDOT  »  AHP  «  f^GA  «  !JCOS(OMEGA«T) 

LOADD  -  2PSL0H«EPD0T 

DLTAEP  -  D?JBS(EPGLOK)  -  SPSLHO 

nr  { DLTAEP ,  XT .  0 . 0 )  Tiffil* 

BETA-(E*ALPA)*(  EPSLOM  -  Z(l)/E  + 

>  PT«ER7(2.*(l,+VKBWIK)/3.»A*EPSLOH)’*UHIT{-LO>J)i>)  ) 
I7(HH.EQ.i.0)THEH 

2DCT(1)  »  E»(  EPDOT  -  2 . * ( l.+VHEMIN) /3. <>DABS(EPK)T) • 

>  (Z(l)-aETA)/y  ) 

ELSE 

7DOT(l)  »  E*{  EPDOT  -  2 .*  ( 1.+VHEK3J) /3.*DA'lS(EPDOT)* 

>  (DABS{2J1)-BETA)/Y>**<HH-1)«(Z(1)-BETA)/Y  ) 

HHDIP 

SLSE 

^•^^(i)  -  e®ep<x;t 
EHDir 
P^TORE 

end 

c 

End  of  UNIAXIAL. FOR 


EXA]^^PLE  OF  AIs  INPUT  FILE: 


SHA  Hysteretic  Model  with  Elasticity  Outside  eO 
855  pai 
28.5x106  psi 
3 

.818 
28.5E06 
855. 

0.818 
3 


0.35 
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APPENDIX  C 


FORTRAN  MAIN  PROGRAM  (al!  subroutines  and  function  subprograms  that  arc  not 
specifically  1  to  the  calculations  being  made  here  arc  given  in  Appendix  H): 


Kaxai  SHEAR. POR 

“«“>  This  pvogr^in  will  carry  out  a  Runga-Kutta  integration  on  the 
atross-straln  equations  of  the  SHA  hysterotic  KOdel  for  SHEAR 
loading.  Hoto  that  inelastic  tehavlor  is  supraeced  beyond  GAMMAO. 
Cclculated  rasults  for  the  stress-atr&tn  response  are  written 
to  user  defined  output  data  files  where  the  data  is  arranged 
in  two  coliuos  with  strain  in  the  first  column  and  stress  in  the 
second.  The  loading  sequence  (hyotereoia)  starts  with  the  first 
line-  of  the  file  and  proceeds  in  order  to  the  laat  line. 

The  resulting  data  files  can  be  used  to  plot  the  curves  via 
GRAPUBR,  which  reads  ASCII  data  file  with  data 
arranged  in  colum  pairs. 


HOMEHCIATOREj 


REAL  COKSTAFTS: 

A  Material  constant  controlling  shape  of  hysteresis 

ALPA  “Ey/IE-Py)  whe'-e  Ey  ic  the  inelastic  Dodulus 
AMP  Anplitude  of  stroll,  input 

E  young ‘ s  modulus 

G  Shear  modulus 

GAMHAO  Strain  beyond  which  inelastic  growth  is  supressed 
FREQ  Frequency  of  cyclic  strain  input 
FT  Material  constant  controlling  size  of  hysterenlB 

N  Overstreas  power  (controls  sharpness  of  transition  to  insl.) 

VHEif  Poisson's  ratio 

y  A^ciai  Stress  where  damping  mechanisms  are  activated 

VS  Shear  Stress  where  damping  mechanisms  are  activated 

IHTSGEP.  CONSTANTS! 

NCyci.E  Mo.  of  cycles  of  oscillation 

NPPCVC  Hr.  o*  points  per  cycle  to  be  used  in  integration 
CHARACTER  STRINGS! 


FIMNAHK 

string  for  filenasoe  assignment  to  a  FORTRAN  unit 

number 

SAL? 

string  for  material  constant  ALFA 

SA 

string  for  materiel  censtent  A 

SG 

string  for  the  shear  modulss  G 

SCO 

string  for  iiaiting  strain  GAMM-AO 

SGP 

string  for  peak  strain  AM? 

SFT 

string  for  material  conotsat  FT 

SK 

string  for  overstrass  power  H 

SVS 

string  for  shear  stress  where  dfeaping  mechanisms 

act ivate : 

TITLE 

descriptive  title  for  run 

VARI.A3IES! 

ENGABS  Energy  absorbed  per  cycle  of  oscillation 

EHGSTO  Energy  stored  (  =  -5  *  icax.  strain  ♦  suxx.  stress) 

ETA  1x380  factor 

SDC  Specific  damping  capacity 


vARiAsr^  APjy.ys! 
GAHMA(Sj 
TAU(K) 
t(K) 


Shear  strain  (index  K  rep.  time) 

Shear  stress  (index  K  rep,  time) 

Shear  stress  passed  fron  INTFUR  (index  JC  rep.  time) 


saSROUTIKESt 

DSRIV 

IKTPDK 

MAXIM 

3TRLES 


Contains  tho  differentiol  eqs.  (invoked  by  IBTPUN) 
Integration  routine  ( «th  order  forward  Runge-Kutta • 
Dstesrmir.es  the  ain.  £  itax.  values  of  an  array 
Counts  the  number  of  characters  in  a  string 


FUNCTION  SUBPPCS.RAMS: 

ERF(X)  Error  function  of  S 

FACT(K)  K  factorial 

SGN(V)  Signum  function  of  V 

UHITjS)  Unit  step  function  of  X 
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c 

c 

IMPLICIT  REflL*8  (A-H,0-Z) 

INTEGER  KBQ 

REAL®8  DERIV,FLOAT,T,TEKD,TOL,Z(20),N 
DIMBKSIOK  GAMMA(lOOl)  ,TAU(1001) 

CHARACTEa*20  PILESAHE 

CHaRACTER®70  TITLE, SYS, SG,SH, SALE, SGO,SGP,SA,SPT 
CHARACTER® ( * )  SARI , SAR2 , SAR3 , SAR4 , SAR5 , SARo , SAR7 , SARP 
PARAMETER  (SARl-'gO»') 

PARAMETER  (SARE-'Ya-') 

PARAMETER  (  SAR3-  ‘  6- '  ) 

PARAMETER  (SAR4-'n-') 

PARAMETER  ( SARS- * alf a- ' ) 

PARAMETER  ( SAR6" '  a- '  ) 

PARAMETER  (SAR7»’fT-') 

PARAMETER  ( SARP- ' gP~ ' } 

COMMON/BLOK1/WC>,OMEGA,YS,E,G,ALPA,K,GAMMAO,A,FT,RAD3 

C 

KEQ  -  1 
C 

T  »  Q.O 
Z(l)  =>  0.0 
C 

C— — >  Interactive  input  of  filenaaea  for  oaterial  data  input  file,  and 
C  for  output  files  contaning  data  and  text. 

C 

WRITE (*,*) ‘Enter  the  naae  of  your  input  data  file.’ 

IFL-9 

1  READ (*, 2) FILENAME 

2  FORMAT(A) 

WRITE  (*,  3 )  FILENAME 

3  FORMAT{  ■  •  ,3X,A20) 

IF(IFL.EQ.9)THEK 

OPEH(IFL, FILE-FILENAME, STATOS-‘UNKNOHN‘  ) 

ELSE 

OPEN(IFL,FILE-FIL2NAMS,STATUS-‘UNKNOWN‘  ) 

ENDIF 

IF  (IFL.EQ,9)  THEN 

WRITE(*.*) 'Enter  naa^  for  the  plot  file  for  TAU  vs.  GAMMA* 

IFL-10 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.IO)  THEN 

WRITE (»,*)' Enter  name  for  the  sunnary  text  file’ 

IFL-11 
GO  TO  1 
ENDIF 
REWIND  09 
REWIND  10 
REWIND  11 
C 

Read  input  quantities  frcn  input  file. 

C 

RFJa)(9,5)TI'nE 
READ(S,  •  (A)  ■  )SYS 
READ(9, ’ (A) ' )SG 
READ(9, ■ (A) ■ )SN 
READ(9, ■ (A) ' jSALF 
5  FORMAT(A) 

C 

C— — «>  Read  in  uniaxial  properties;  convert  to  shear  properties  later. 

C 

RBAD(9,*)2,VHEW 

READ(9,*)y 

READ(9,»)ALFA 

READ(e,*)N 

C 

C— “>  Interactive  input  of  other  Material  paraasters,  strain  eunplitude, 
C  and  number  of  cycles  of  loading  tc  be  used  in  calculations. 

C 

WRITE (*,*)■  Enter  GAMHAO ’ 

READ(S,*)6AMHA0 

WRITE Enter  character  string  for  GAMHAO' 

READ(5,  •  (A)  ■  )SG0 
WRITE(*,*)'  Enter  A  AND  FT' 

READ (5,*) A, FT 

WRITE (*,*)'  Enter  character  string  for  a ' 

READ(5, ' (A) • )SA 

WRITE (»,*)'  Enter  character  string  for  fT' 

READ(5, ■ (A) ’ )SFT 
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'  Eutst  FREO  oxid  AMP  Of  giM&7  strain  loading' 
RSAD(5,»)PR2Q,i!HP 

WRITE  (*,*)'  Enter  cUarac-tar  string  for  peak  shear  strain  ’ 

READ{5, • (A) • )SGP 

WRITE  ( * ,  ♦ )  ■  Enter  Kunhar  of  Ch'cleo  and  l.'ueber  of  Points  oer  t'ycle ' 
READ  ( 5 ,  • )  NCYCte  ,  ar-PCEC 
c 

Dsteraine  length  of  '.ario’io  character  strings  for  iater  use. 

C 

CALL  STR3J3H{SGO.T;3MO,IESeO) 

CALL  STRLKK(SG?,IBS»3I'.IiESGFj 
CALL  STltL5S(SVS,IBSXP,IESyS) 

CALL  STaLEH(RG,IBSG,IESG) 

CALL  SSRLEH(SH,IDSH,IESH) 

CALL  STRLEHfSALF.lBSALF.KSAL?) 

CALL  S'rRLEHi5  ',1BSA,IESA) 

CALL  STRLEB(SFT,I3SFT,I2SPT) 

CALL  STRLEH(SAR1,IBSAR1,I2SAR1) 

CALL  STRLEH(SARP,IBSARP,IESASP) 

CALL  STRLBK(SAR2,1BSAR2,IESAR2> 

CALL  STRLEH(SAR3,1BSAR3,IESAR3) 

CALL  STRLBH(SAR4,IBSAR4;IESAR4) 

CALL  STRLEH(SAP^,IBSARS,IBSARS) 

CALL  STRLE»(SAR6,IBSAR6,IESAR4) 

CALL  STRLEM(SAR'’,1BSAR7, J2SAR7) 

n 

ISDMl  -  lESARi  4  lESGO  +  2  +  lESARP  +  lESGP 

ISUKl  -  IESAR2  +  lESY  +  2  +  IESAR3  +  lESE 

ISCM3  «•  tbsJiR4  ^  ibsH  +  2  +  IESAR5  +  lESALF 

ISOMi  "  IESAR6  +  lESA  +  2  +  IESAR7  +  lESFT 

C 

C“«“>  Calculate  the  quantity  Pi"3.141592S7. . . .  and  other  paraaeters. 

C 

PI"DACOS(-1.0D00) 

RAD3“3.*».5 
G  «  2/{2-*(l+VlJEW)  ) 

YS  -  Y/HAD3 
OMEGA  »  2.*PI»FREQ 
PRINT  PI-',PI,  '  Q«EGA»' .OMEGA 
KSTEPS  -  KCYCLE»NPPCYC 
PRINT  KSTEPS*' , KSTEPS 

PERIOD  >  1./PR2Q 
DELT  -  PERIOD/NPPCYC 
PRINT  DELT- '.DELT 

C 

C“»=*->  Set  initial  conditions 
C 

TAD{1)  ~  0. 

GAMKA(l)  -  0. 

C 

WRITE  (*.99999) 

99939  FORMAT ( 4 X, 'ISTEF'.SX. 'TIKE'.9X, 'Zl'.llX, 'Z2•^ 

C 


Carry  oat  the  numerical  integratior 


C 

C 


C 


DO  10  ISTEP  ~  1-NSrZPS 

CALL  INTFUN{Z,T.DELT.NEQ) 

TAD(ISTEP+1)  »  S;l) 

GAHHA(I3TEP+1)  -  AMP  •  DSIN(OfuiGA*T) 

WRITE  (*,  '  (SX,I6,4(3X.E10.4) )  '  )ISTEP  GAHMAJISKP+l )  ,S  (  H 
10  CONTINUE 


50 


CALL  KASIM (TAU.HSTEPS.TADHIN.TAUMAX) 

ENGABS  »  0. 

ISTART  -  HPPCyC/4,+1 
lEND  =  ISTART  +  NPrCYC-l 
DO  50  I-ISTART. lEHD 

ENGABS  »  EKGMS  +  .5  »  (T.\D(I+1)  +  TAD(I)j  « 

■  (GAMMA(l-ii)  -  GAMMa(I)) 

CONTIKBB 

EKGSTO  «  .5*T7.UMAX*AH? 

ETA  =  ENGABS/ ( 2. *PI*EHGSTO) 

SDC  ETA*2*P2 


c*—- >  Write  the  results  to  output  aata  file  (unit  10)  and  to  output 
C  text  file  (unit  11). 

C 


WRITE(10. ' flX.2(E16.8.  ' ' ) ) ' ) (GAHKA(I) ,TAU;i) . I-1.K3TEPS+1 ) 

WRITS(H. '  (IX.  "0  ",I2."  *■ '.A.A,2X.A.A . )') 

>  ISOf41,SAai,aG0(IBSG0iISSG01,SABP.SGP(I3SGPlI?.rCP) 


DTRC-SME-91/34 


WRITS(il,  •  (IX,  •  "  ,12, . ,&,&,2X,A,£, . )  ’ ) 

>  ISUH2 , SAR2 , SYS ( IBSYS  t lES YS ) , SAR3 , SG ( IBSG  S lESG ) 

RRITEdl,  ■  (IX,  ■  '2  •  ',12, . ,A,A,2X,A,A, . )  ’  ) 

>  ISnM3,SAR4,SN(IBSHlIESN),SARS,SALP(IBSALrjIESALF) 

VfRITEdl,  •  (IX,  •  ‘3  •  •  ,12, . ,A,A,2X,A,A, . )  '  ) 

>  ISnM4  ,  SAR6 ,  SA(  IBSAl  lESA)  ,  SAR7  ,  SrY(  IBSFT I  lESPT) 

WRITE(li, ‘{IX,  • '4  28  "ETA- '  ’  ,E9 .4,2X,  ’ ‘SDC-' ‘  ,E9 . 4, . )')ETA,SDC 

PRINT  LOSS  FACTOR-' , ETA,  •  SDC-',SDC 

C 

1000  FORMAT (A70) 

1001  FORMAT(2A60) 

1200FORMAT('  ' , 3X, ‘ 1 ’ ,5X, 17 ) 

1201  FOHHATC  •,3X,  •2',SX,I7) 

1500  FORKA'P(*  •,5X,E10.4,SX,E10.4) 

1501  FOHMATC  '  ,52,E10.4,5X,E10.4,5X,S10.4) 

STOP 

END 


SOBROOTINE  DEaTV(Z,T,ZDOT) 

IMPLICIT  REAL'S  (A-H,0-2) 

INTEGER  NEQ 

REAL'-e  OKEGA,T,i;(20),SOOT(20),NH,LOAD» 

COMHCM/BLOK1/AMP,OK'GA,YS.-S,G,ALFA,HH,GAMHAO,A,FT,RAD3 

C 

GAMMA  <»  AMP  *  OSIN(CMEGA»T) 

GEMDCT  -  AMP  '  OMEGA  *  DCOS(C»a:GA*T) 

DLTAGA  -  DABS  (GAMMA)  •  GAMKAC 
IF ( DLTAGA . LT . 0 . 0 ) THEN 

BETA»(1./3.)*(E*AL?A)*(  GAMMA  -  Z(l)/G  + 

>  RAD3*PT»ERF(A*GAhMA/RAD3)*imiT( -GAKMA'GAMBOT)  ) 
lP(KN,EQ.l.rt)T8ES 

ZDOT(l)  -  G>'{  GRMDOT  -  DABS(GAHBOT;»{2(l)-BETA)/YS  ) 
ELSE 

2DCT(1)  -  G»(  GAMDOX  -  PABS(GAKDOT)* 

>  (DABS(2(l;-fiE'fi.)/YS)’-*(KN-l)  »  (Zd)-BBTA)/YS  ) 
ENDIF 

ELSE 

SDOT(l)  G'GAMDOT 
BN0i.F 
RL-niRN 
END 

C— SHEAR.  FOR 


EX  AMPLEOFAN  INPUT  HLE: 


SHA  Hysteretifs  Model  for  Shear  with  Elasticity  Outside  gO 
494  psi 
10.56x156  psi 
3 

.818 

28.5E06  0.35 

855. 

0.8x6 

3 
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APPENDIX  D 


FORTRAN  MAIN  PROGRAM  (all  subroutines  and  function  subprograms  that  are  not 
specifically  related  to  the  calculations  being  made  here  are  given  in  Appendix  H): 


Pfogram  Hases  SDUA.FOR  (for  £tialn  dependent  Qni^ial  response) 

““>  ‘fhis  progran  is  different  froB  DNIAXtAL.FOR  in  that  it 

repeats  the  cyclic  strain  application  over  a  specified  range 
for  a  given  nucber  of  peak  strains,  and  it  calculates  the  loss 
factor  associated  with  each.  Rungo-Kutts  integration  will  be  applied 
to  tno  stress-strain  eguaticns  of  the  axial  SKA  hyeteretic  oodel. 
Hote  that  ine'.astic  behavior  is  supressed  beyond  ePSLHO. 

Files  will  be  >ensrated  which  contain  data  for  the  plotting 
of  Peak  Strain  vs.  Energy  Abrorbed,  Peak  Strain  vs.  Loss  Factor, 

Peak  Eguivalent  Srain  vs.  Energy  Abt orbed,  and  Peak  Equivalent 
Strain  vs.  Loss  Factor.  Pile  Pomat  fill  be  that  of  GRAPHER; 

ASCII  data  is  arranged  in  coluxm  pairs.  Data  is  arranged  in 
two  vertical  colunns  Nith  strain  (or  equivalent  strain)  being 
in  the  first  coluxm  and  lose  factor  (or  energy  absorbed)  being 
in  the  second  coiuiui.  A  short  set  of  date  sunxuirizing  the  run 
is  written  to  a  user  defined  text  file. 


K -KM  liAWiAM  M 


ROMEHCLATUREs 


P£AL  CCNSTAKTS; 

A  Haterial  constant  controlling  shape  of  hysteresis 

ALFA  -Ey/(E-Ey)  where  Sy  is  the  i.^elastic  eodulus 
AMP  Amplitude  of  strain  input 

E  young ' s  zodulus 

SPSIRO  Strain  beyond  which  inelastic  growth  is  supressed 

EPSPXl  Kininuxi  peak  strain 

EPSPK2  Haxiaiur  peak  strain 

FREQ  Frequejicy  of  cyclic  strain  input 

FT  Material  constant  controlling  size  of  hysteresis 

N  Oveistress  power  (controls  sharpness  of  transition  to  inel.) 

VNEW  Poisson's  ratio 

VHEW7R  Inelatic  Poissen  coefficient 

Y  Stress  where  daaping  isechanisDs  are  activated 

INTEGER  CONSTANTS! 

NCYCLE  No.  of  cycles  of  oscillation 

NPPCyc  No.  of  points  per  cycle  to  be  used  in  integration 
NINC  No.  of  in^renents  for  the  ra^ge  of  peak  strain 

CHARACTER  STRINGS: 

FILENAME  String  for  filenanc  assignment  to  a  FORTRAN  unit  number 

SALF  string  for  material  constant  ALFA 

SA  string  for  material  constant  A 

SE  string  for  Young's  modulus  E 

SEO  string  for  limiting  strain  EPSLNO 

SFT  string  for  material  constant  FT 

SN  string  for  cverstress  power  H 

SY  string  for  stress  where  damping  xiechanisms  activate:  Y 

TITLE  descriptive  title  for  run 


VARIABLE  ARRAYS! 
EHRGAB(I) 
ENRGST{I) 
EPSLON(X) 
2TA(I) 
PEAKST(I) 
PKEQST(I) 
STRESS (K) 

Z(K) 


energy  absorbed  in  ore  cycle  (1  recresents  position) 
Energy  absorbed  in  one  cycle  (I  represents  positron) 
Axial  strain  (index  K  rep.  time) 

Loss  factor  (I  represents  position) 

Peak  strain  (index  I  represents  position) 

Peak  equivalent  etrain  (index  I  represents  position) 
Axial  stress  (index  X  rep.  time) 

Axial  stress  passed  from  INTPUN  (index  K  rep.  trme) 


SUBROUTINES ! 
DERTV 
IKTPUN 
MAXIM 
STRLEN 


Contains  the  differential  eqs.  (invoked  by  IHTFJH) 
Integration  routine  («th  order  forward  Hunge-Kutta) 
Determines  the  min.  i  max.  values  of  an  array 
Counts  the  nutiber  of  characters  in  a  string 


FUNCTION  SUBPROGRAMS! 

BRP(X)  Error  function  of  X 

FACT(K)  K  factorial 
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SGN(V)  SignuB  function  of  V 

nNIT(X)  Unit  step  function  of  X 


IMPLICIT  REAL*8  (A-H,0-S) 

INTEGER  HBQ 

REAL*8  DERIV,rLO.\T,T,TEND.Z(20),H 
DIMENSION  EPSLON(1001),STRESS(1001), 

>  PEAXST( 200 ) ,ETA( 20C) ,EHRGAB( 200 ) ,  ETOGSTj  200 ) ,PKEQST( 200 ) 

CHARACTER* 20  FILENAME 

CHARACTER*70  TITIi ,  S Y ,  SB ,  SN ,  SALP  ,  SEO  ,  SA,  SPT 
CaARACTER*(*)  SARl,SAa2,SAP.3,SAR4,SAR5,SAH6,SAR8 
PARAMETER  ( SARI- ' e0» ' > 

PARAMETER  (SAR2-‘Y-‘ > 

PARAMETER  {SAR3-’E-' ) 

PARAMETER  (SAR4-‘n-'i 
PARAMETER  (SARS-'n-*) 

PARAMETER  ( SAR6- ‘ a- ‘ ) 

PARAMETER  (SAR8~'fT-'» 

COMHCN/DLOKl/AHP, OMEGA, Y,E,VNEW,VNEKIH,EPY, ALFA, N.EPSLNO,  A,  PT 
NEQ  -  1 

Interactive  input  of  filenaises  for  the  luitericl  data  input  file 
and  for  output  files. 

WRITE  (*,*) 'Enter  the  nasks  of  vour  input  date  file.' 

IFL-9 

1  READ (*, 2) FILENAME 

2  FORMAT(A) 

WRITE  (*,  3 1  filename 

3  FORMATl'  •,3X,A20) 

I?(I.»L,Ee.9)THEH 

OPEN  ( IFL,  FILE-FILENAME ,  STATUS- '  UNKNOWN '  ) 

ELSE 

OPEN  (IFL,  FILE-FILENAME,  STATUS- 'UNKNOWN'  ) 

XNDIF 

IF  (IPL,EQ.9)  THEN 

WRITE (*,*)' Enter  nane  of  plot  file  for  Loss  Fac.  vs.  Peak  Str. ' 
IFL-IO 
GO  TO  1 
ENDIF 

IP  (IPL.EQ.IO)  THEN 

WRITE (*,») 'Enter  naase  of  the  GRAPHER  legend  text  file' 

IFL-ll 
GO  TO  1 
ENDIF 

IP  (IFL.EQ.ll)  THEN 

KRITE(*  ,*) 'Enter  filenasie  for  Energy  Absorlksd  vs.  Peak  strain.' 
IFL-12 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.12)  THEN 

WMTE (*,*)' Enter  filenane  for  Loss  Factor  vs.  Peak  Eqiv.  strain' 
IFL-13 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.13;  THEN 

WRITE (♦,»)' Enter  filenaxte  for  Energy  Abe.  vs.  Peak  Eqiv.  Strain' 
IFL-14 
GO  TO  1 
ENDIF 
REWIND  09 
REWIND  1C 
REWIND  11 
REWIND  12 
REWIND  13 
REWIND  14 

■— »>  Reed  input  quantities  frcn  input  file, 

READ  (?,  3)  TITLE 
RPAD{9, ' (A) ' )SY 
REAa(9, ' (A) ' )dE 
READ(9, ' (A) ' )SN 
READ(9, '(A)')SALF 
5  FORHAT(S) 

READ(9,*)E,VHEW 
READ(9,*)Y 
READ ( 9, *J ALFA 
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HBAD(9,*)H 


Interactlvt  itput  vt  other  saaterial  para».jter8,  ottair  oispljtuda, 

C  JUKI  nunoar  off  cyel«s  of  leading  to  t>a  used  in  calcs.  \ati£  ">8 

C 

WRISE {*)*)*  Enter  EPSiao ' 

REAJ>{5,»)ErCUI0 

»mZTE  <*.*)'  EPSLNO-  ■  ,  EP.9LH0 

WslTE(*j»)‘  Eater  chorj-oter  *tring  for  EPSrso* 

READ(5,‘(;a)')SE0 
WRiyB{*,’(lZ,  ••  SE0-'\A}',SE0 
»31ITE(*,*}‘  Enter  A  A»K>  FT* 

REAC{5,*)A,?T 

WftITa(*,*)’  A''’, A,'  PT-*,r^ 

WRITE}*,*; •  Eater  character  string  for  a* 

REAE(5,'(A)‘)<5A 

'WRITE}*, ‘{iX,  •  '  SA-' '.A) ‘jSA 

SPJl'E } * , * )  •  Enter  character  t-triag  for  fT‘ 

8JsAE(5,’}A)‘)SPT 

WRITE}*,  •}1X,'‘  SPT-J' '.AJ  •  )f7T 

WRITE}*,*) •  Enter  PREQ,  EPSPXl-  S  EPSPK7  }ain  and  anx  peak  atr. ) ' 

READ } 5 , • ) PHEO, EPSPKl , EPSPKi 

WRITE;*,*)’  ?HE0-’,FR2Q, •  EPSPKl-’ ,BP6PK1, ’  EPSPKS-' ,SPSPK2 
WRITE}*,*)*  Enter  Ko,  of  CyolAp  and  Ko.  of  Feints  per  Oycl'i' 

P,2?i: }  5  ,  * )  HCi'CLE ,  WPPCYC 

WRITE  }*,*)•  HC YCliE-  ’ ,  KCYCXE  ,  '  NPPCyc«  •  ,  tJopCYC 

WRITE}*,*)'  Enter  NIHC  JHo.  of  incresents  bat.  peak  strains)' 

H2AD;5,*)BtKC 
WRITS}*,*)'  SIKC-'jHIBC 
C 

Dotaraiae  length  ot  'various  character  strlnos  for  later  I'se. 

C 

CALL  STaL2H(SP0,iB3E0,I8SE0) 

CALL  arRLEK}3Y,IBSY,IESY) 

CALL  STRL?.K;SE,IBSE,iESE) 

CALL  STRLEB}SN,IBSK,IESH) 

CALL  STKLEB}SVJr,IBSALP,IESALP} 

CALL  S?RLEH}SA,IBSa,iESA) 
a\LL  S'm£H}SPT,IBSPT.2ES7T) 

CALI,  STRLBH }  SARI ,  IBSARl ,  XESARi ) 

CALL  8'JRLEV}SAR2,IBSAR2,IESAR2) 

CALI,  STHLEa}SAR3,IBSAR3,IESAR3) 

CALL  S1RLEH(SAR4,IBSAR4,XSSAR4> 

CALL  STRLEH }  SAR5 ,  IBS.AR5 ,  lESARh ) 

CALL  STRLEH (SAR6,1.!3SARS,IB<?AR6> 

CAIL  STRLEH }SAS8,xB2ARe,IESAR0) 

C 

XSUMl  •>  lESABi  lEStO 

1SUM2  «  1ESAP.2  +  lESY  +2  -vlESARS  +  T2SE 

IStfM3  ■*  IHSAR4  +  ISSH  +  2  +  IESAR5  +  lESAL? 

ISUH4  “  lESARC  +  ISSA  +  2  +  IESAE8  +IESFT 

C 

Calculate  the  quantity  Pi-3. 141592(17. ,  and  other  pjsranetero 

PI-D.ACOS ;  -1 ,  OBOC ) 

Of»EGA  -  2.*?1*FP.EC' 

PRIHT  *,'  PI-', PI.'  <»3GA»=' .(SEGA 
VME'WIH  -  .5  -  AI.FA*}.5-VKEW)/}J.+ALPA) 

EPy  -  Y/a 

PRIin-  *,'  VHEWIN-',VHEWT(I,  •  EPY»’,EPY 
NSTEP3  -  HCYCniiE'HPPCYC 
NSTPri  -  HSTEPS  r  1 
PRIKT  *,  '  NSTSPS-' , HSTEPS 
PERIOD  “  l./?REQ 
DELT  *■  PERIOD/SPPCYC 
PRIHT  *,'  DELT-',DBLT 
EPSIKC  -  }S?S?K2-EPSF51)/NI.HC 
ISTART  -  fIPPCyC/4.  +  1 
lEND  -  ISTART  +  HPPCfC-l 
C 

C'«-— >  Set  initial  conditions. 

C 

STRESStl)  -  0. 

EPSLON}!)  -  0. 

C 

C— — >  Pet  up  the  loop  for  carrying  out  integration  for  each  step  ir.  peak  strain. 


KIRCPl  -  KINC  +  1 
DO  7777  J-1,HIHCP1 

■>  Clear  the  arrays  used  in  the  loop  narked  by  30  at  the  continue  stateaent 
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Reinitialize  T  and  Z ( 1 ) . 

DO  10  K-1,NSTPP1 
EPSLON(X)-0.0 
10  STRESS(K)-0.0 

T  -  0.0 
Z(l)  -  0.0 

Update  the  peidc  strain 

AMP  -  EPSPKl  +  EPSXMC»(J-1) 

PEAKST(J)  -  AHP 
PKEQST(J)  -  (1+VHEW)*AMP 

“»">  Ceurry  out  the  integration  for  the  current  peak  strain 

DO  30  ISTEP  -  l.HSTEPS 

CALL  INTFUH(Z,T,DELT,NEQ) 

STRESS (ISTEPtl)  ••  Z(l) 

EPSLON(ISTEP+l)  -  AHP  *  DSIK(OMEGA*T) 

30  CONTIHUE 

“•«““>  Coi^>ute  danplng  for  the  current  peak  strain. 

STRMIH  -  0. 

STRHAX  -  0. 

CALL  HAXXH( STRESS, NSTPPl, STRMIH, STRHAX) 

EHGABS  -  0. 

DO  50  I-ISTART,IEND 

EHGABS  -  EHGABS  +  .5  •  (STRESS(I+1)  +  STRESS(I))  * 

>  (EPSLOH(I+l>  -  EPSLOH(I)) 

50  COHTIHUE 

EHGSTO  -  .5*STRHAX«AHP 
EHRGAB(J)  -  EHGABS 
EHRGST(J)  •  EHGSTO 
ETA(J)  -  EHGABS/ ( 2. *PI*EHGSTO) 

WRIT2(», ' {3X,I3,3{3X,E10.4)) •)J,PEAKST(J),ENRGAB(J),ETA(J) 
7777  COHTIHUE 

CALL  HAXIH(ETA,HIHCP1,ETAHIH,ETAHAX) 


Write  renults  to  output  data  file  (unit  10,  12,  13,  and  14)  and  to 
output  text  file  (unit  11). 


WRITE(10, •(1X,2(E16.8, *',’•)) ' ) (PEAKST(I) ,ETA(I) ,I-1,HIHCP1) 
WRITE(11,  •  (IX,  "0  ",12, . ,A,A, . )’) 

>  ISUH1,SAR1,SE0(IBSE0IXESE0) 

WRITE(11,  ■  (IX,  "1  ",I2, . ,A,A,2X,A,A, . )  ‘  ) 

>  ISUH7 , SAR2 , SY ( IBS Y : lES Y ) , SAR3 , SE ( IBSE : lESE ) 

WRITE(11,  '  (IX,  •  '2  "  ,12, . ,A,A,2X,A,A, . )  '  ) 

>  ISUH3  ,  SAR4  ,  SH  ( IBSH  :  lESH )  ,  SAR5  ,  SALT  ( IBSALE :  lESALF ) 

HRITE(11,  •  (IX,  "3  ",12, . ,A,A,2X,A,A, . )  '  ) 

>  ISUH4,SAR6,SA(IBSA:IESA), 

>  SAR8,SFT(IBSFT:IESFT) 

WRITE(11,  '  (IX,  "4  38  "ETAjoax-"  ,E9.4,2X,  "ETAain-"  ,E9 .4 , . )') 

>£TAHAX,STAHIH 

WRITE  (12,  •  (1X,2(E16.8,  ",  "  )  )  '  )  (PEAKST(I)  ,EHRGAB(I) ,  I-l  .HIHCPl ) 
WRITE(13, • (1X,2(E16.8,  ",  "  ) ) ' ) (PKEQST(I) ,ETA(I) , I-1,HIHCP1 ) 
WRITE(14,  '  (1X,2(E16.8,  ","))')  (PKEQST(I)  ,EHRG7iB(I) ,  I-l  ,HIHCP1 ) 
PRIHT  •  ,  ■  ETAHAX- ' ,  ETAH70C ,  '  ETAHIH- '  ,  ETAHIH 


1000  FORHAT(A70) 

1001  FORHAT(2AeO) 


1200  FORHAT( ' 

1201  FORHAT( ' 

1500  FORHAT( ' 

1501  FORHATC 
STOP 
EHD 


■,3X,'1',5X,I7) 

■,3X, •2',5X,I7) 

' ,SX,E10.4,SX,E10.4) 
',5X,E10.4,5X,E10.4,5X,E10.4 


) 


SOBROBTIHE  DERIV(Z,T,ZDOT) 

IHPLICIT  REAL*8  (A-HfO-Z) 

R£AL*8  OMEGA,T,Z(20),ZDOT(20),HH,LOADD 

COHHOH/BLOKl /AHP , OMEGA, Y , E , VNEW, VHEWIH, EPY, ALFA , HH , EPSLHO , A , FT 

EPSLOH  -  AHP  *  DSIH(OHEGA»T) 

EPDOT  -  AHP  •  OMEGA  •  DCOS(OHEGA*T) 

LOADD  -  EPSLOH  *  EPDOT 
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DMAEP  -  DABS(EPSI.CN)  -  EPSUJO 
C 

IP ( DLTAEP • LT . 0 . 0 ) THEN 

IF(DABS(EPSLON) .GE.EPY)TH£N 

DELTA  -  .5  -  DABS(Z(l)/EPSLON)«(.5-VNEW)/E 
ELSE 

DELTA  -  VNEW 
EHDIP 

BETA-(E»ALPA)*(  EPSLON  -  Z(l)/E  + 

>  PT»ERF{ (2.*(l.+DELTA)/3. )*A*EPSL0N)*UHIT(-L0ADD)  ) 
IP(NH.EQ.1.0)THEH 

ZDOT(l)  -  E*(  EPDOT  -  (2.*(i.+VKEWI»)/3.  )*t)ABS(EPDOT)* 

>  (Z(l)-BETA)/y  ) 

ELSE 

ZDOT(l)  -  E*(  EPDOT  -  (2.*{l,+VNEWIN)/3.  )*I>ABS(EPDOT)* 

>  (DABS(Z(1)-BETA)/Y)**(HN-1)*(Z(1)-BETA)/Y  ) 

ENDIP 

ELSE 

ZDOT(l)  -  E»EPDOT 
EHDIP 
RETDRK 
END 


EXAMPLE  OF  AN  UT  FILE: 


SMA  Axial  Hysteretic  Model  with  Elasticity  Outside  eO 
855  psi 

28.5x106  psi  Vnew".35 
3 

.818 

28.5E06  0,35 

855. 

0.818 

3 
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APPENDIX  E 


FORTRAN  MAIN  PROGRAM  (all  subroutines  and  function  subprograms  that  are  not 
specifically  related  to  the  calculations  being  made  here  are  given  in  Appendix  H): 


>.««>_«>  progran  None:  SDSH.POR  (for  Strain  Bependent  Sflear  response) 

■“>  This  progran  is  different  frora  UNIAXIAL. FOR  in  that  it 

repeats  the  cyclic  strain  application  over  a  specified  range 
for  a  given  nunber  of  peak  strains,  and  it  calculates  the  loss 
factor  associated  with  each.  Runge-Kutta  integration  will  be  applied 
to  the  stress-strain  equations  of  the  SKA  hysteretic  nodel. 

Note  that  inelastic  behavior  is  supressed  beyond  GAHHAO. 

The  applied  loading  is  that  of  a  pure  state  of  shear  stress. 

He  specify  a  neccuiujD  peak  strain  to  be  used  in  the  sinusiodal 
application  of  shear  strain.  The  resulting  cur/es  for  stress  and 
strain  will  be  stored  in  data  arrays  and  then  will  be 
used  to  calculate  the  loss  factor  (danping).  The  data  file  will  also 
be  in  a  fcmat  conpatible  for  plotting  with  GRAPHRR. 

ASCII  data  is  arranged  in  coluim  pairs.  Data  is  arranged  in 
two  vertical  coluxins  with  strain  (or  equivalent  strain)  being 
in  the  first  colunn  and  loss  factor  (or  energy  absorbed)  being 
in  the  second  coluim.  A  short  set  of  data  suonarlzing  the  run 
is  written  to  a  user  defined  text  file. 


NOMENCLATURE: 

REAL  CONSTANTS: 

A  Material  constant  controlling  shape  of  hysteresis 

ALFA  -Ey/(E-Ey)  where  Ey  is  the  axial  inelastic  nodulus 
AMP  Anplitude  of  strain  input 
E  Young's  nodulus 

G  Shear  modulus  {-E/(2( 1+VNEW) ]} 

GAHHAO  Strain  beyond  which  inelastic  growth  is  supressed 

GAHPKl  Minimum  peak  strain 

6AMPK2  Maximum  peak  strain 

FREQ  Frequency  of  cyclic  strain  input 

FT  Material  constant  controlling  size  of  hysteresis 

N  Overstress  power  (controls  sharpness  of  transition  to  inel.) 

VNEH  Poisson's  ratio 

Y  Axial  stress  where  damping  mechanisms  are  activated 

YS  Shear  stress  where  damping  mechanisms  are  activated 

INTEGER  CONSTANTS: 

NCYCLE  No.  of  cycles  of  oscillation 

NPPCYC  No.  of  points  per  cycle  to  be  used  in  integration 
NINC  Ho.  of  increments  for  the  range  of  peak  strain 

CHARACTER  STRINGS; 

FILENAME  string  for  filename  assignment  to  a  FORTRAN  unit  number 

SALF  string  for  material  constant  ALFA 

SA  string  for  material  constant  A 

SG  string  for  Young's  modulus  G 

SGO  string  for  limiting  strain  GAHMAO 

SFT  string  for  material  constant  FT 

SN  string  for  overstress  power  N 

SYS  string  for  shear  stress  where  damping  mechanisms  activate 
TITLE  descriptive  title  for  run 

VARIABLE  ARRAYS: 

ENRGAB(I) 

ENRG3T(I) 

GAHMA(K) 

ETA(I) 

PEAKST(I) 

PKEQST(I) 

TAU(K) 

Z(K) 

SUBROUTINES: 

DERIV 
INTFUN 
MAXIM 
STRLEN 


Energy  absorbed  in  one  cycle  (I  represents  position) 
Energy  absorbed  in  one  cycle  (I  represents  position) 
Shear  strain  (index  K  rep.  time) 

Loss  factor  (I  represents  position) 

Peak  strain  (index  I  represents  position) 

Peak  equivalent  strain  (index  I  represents  position) 
Shear  stress  (index  K  rep.  time) 

Shear  stress  passed  from  INTFUN  (index  K  rep.  time) 


Contains  the  differential  eqs.  (invoked  by  INTFUN) 
Integration  routine  (4th  order  forward  Runge-Kutta) 
Determines  the  min.  t  max.  values  of  an  array 
Counts  the  number  of  cheuracters  in  a  string 
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FUNCTIOH  SUBPROGRAMS I 
ERF(X)  Error  function  of  X 

PACT(K)  K  factorial 

SGH(V)  SignuB  function  of  V 

UHIT(X)  Unit  atep  function  of  X 


IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  NEQ 

REAL*8  DERIV,FLOAT,T,TEND,TOL,Z{20),H 

DIMENSION  GAHHA( 1001), TAn( 1001 ),PEAKST(201),ETA(201),PKECST(201), 
>  ENRGAB(201) ,ENRGST(201)  . 

CHARACTER*20  FILENAME 

CHABACTER»70  TITLE, SYS,SG,SK,SALP,SGO,SA,SPT 
CHARACTER*(*)  SARI , SAR2 , SAR3 , SAR4 , SAR5 , SAR6 , SAR7 
PARAMETER  ( SARI- '  gO- '  ) 

PARAMETER  ( SAR2- '  YS- '  ) 

PARAMETER  (SAR3-'G«') 

PARAMETER  ( SAR4- ' n- ' } 

PARAMETER  (SAR5-'a-') 

PARAMETER  (SAR6-‘a-') 

PARAMETER  { SAR7- '  f T- '  ) 

CCMMOH/BLOXl /AMP  ,  OMEGA ,  YS  ,  E ,  G  ,  ALFA,  N  ,  GAHMAO  ,  A,  FT ,  RA03 
NEQ  -  1 

■“>  Intcractve  input  of  filenanea  for  the  naterial  data  input  file 
and  for  output  files. 

WRITE (*,*) ‘Enter  the  nane  of  your  input  data  file.' 

IFL-9 

1  READ  {*,2)  FILENAME 

2  FORMAT(A) 

WRITE  C,  3 )  FILENAME 
3PORMAT(‘  ',3X,A20) 

IP{IPL.EQ.9)THEH 

OPEN(IPL,FILE»PILENAME,STA‘rJS-‘ONKNOWH‘ » 

ELSE 

OPEN  ( IFL ,  FILE-PILENA.ME ,  STATUS-  ‘  UNIQJOWH  ‘ ) 

ENDIF 

IP  (IPL.EQ.9)  THEN 

WRITE {♦,•) 'Enter  file  name  for  Loss  Factor  vs.  Peak  Strain' 
IFL-10 
GO  TO  i 
ENDIF 

IF  (IFL.EQ.IO)  THEN 

WRITE {»,*) 'Enter  nane  for  the  suxmary  file' 

IFL-11 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.ll)  THEN 

WRITE (*,*)' Enter  filensse  for  energy  absorbed  vs.  peak  strain' 
IPL-12 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.12)  THEN 

WMTE(*,*) 'Enter  filenane  for  loss  factor  vs  peak  equiv.  strain' 
IFL-13 
GO  TO  1 
ENDIF 

IF  (IPL.EQ.13)  THEN 

WRITE (•,♦) 'Enter  filename  for  energy  abs.  vs  peak  equiv.  strain' 
IFL-14 
GO  TO  1 
ENDIF 
REWIND  09 
REWIND  10 
REWIND  11 
REWIND  12 
REWIND  13 
REWIND  14 

“>  Read  input  quantities  from  input  file. 

READ (9, 5) TITLE 
RBAD(9, '  (A) ' )SYS 
READ(9, ‘ (A) ' )SG 
READ(9, '(A)‘)SN 
READ(9,'(A)')SALF 
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HEA0(9,*-)2,VK2W 

REM>(9,*)ALI'X 

REAn(9,«)H 

C 

C— ><x»:>  intoractive  input  of  other  aatariai  paraaeters,  strain  awplitude, 
C  and  nu!9ber  of  cycles  of  loading  to  be  used  in  calcaletions • 

C 

NRITE(*,*)'  Enter  GXMMAO' 

REA0(5,«)GAHKXC 

MRITE(*,*)'  GAiOOLO-' ,GAMMA0 

WRITE {*,•)*  Enter  character  string  for  GAMMAO* 

RKAI>(5,  •  (A)  ‘  )SGO 

WRITB(*, ' (IX. ' 'SGO-' ’.A) "jSSO 

\«RITS(*,*)'  Enter  A  AHD  P?' 

R2AD(5,»)A,yT 

:^TE(*,*)'  A-*,A,  '  PT“’,PT 

WRITE (*,*)'  Enter  character  string  for  a' 

READ<S,  '  (A)  '  )SA 

«RITE(*, • (IX, ’ ’SA-* ■ ,A5 ' )SA 

WRITE (*,“)'  Enter  character  string  for  f T 

READ(5, ■ (A) ’ )SPT 

KRJTE',*,  •  (IX,  ‘  ’SPT-'  *  ,A)  '  )SPT 

VJRITE(»,»)'  Enter  FRBQ,  GAHPKi,  S  GAHPK2  (Min  £  Max  Peak  Strains)' 
REAS  ( 5  ,  * )  FREQ ,  GAKPKl ,  GAMPK2 

WR1TE(*,*)'  FREQ-' ,PREQ,  ‘  GAMPKl- '  .GAKPAl  ,  '  GAMPK2-' ,GAMPK2 
WRITE {*,«)'  Enter  Hunber  of  Cycles  and  Utuaber  of  Points  per  Cycle' 
READ ( S , • ) HCYCLE , HPFC YC 

WRITE (*,»)'  HCYCLE- ' , HCYCLE , '  HPPCYC- ' , HPPCYC 

5<RITE(*,»i'  Enter  KING  (So,  of  increicents  between  peedc  strains)' 
READ(S,*)SrBC 
HRITE(«,«)'  HINC-‘,HI»C 
C 

C— — >  Detentine  length  of  various  character  strings  for  later  use, 

C 

CALL  STRLBH(SGO,IBSGO,IESGO) 

CALL  STPiEN(SyS.IBSYS,IESYS) 

CALL  STSLEK(SG,IBSG,IESG) 

CALL  STRLEN(SH,I3SS,IESH) 

CALL  SXRLEK(SALF,IBSAL?,IESALP) 

CALL  STRLEN(SA,1BSA,IESA) 

CALL  STRLSN(SPT,IBSPT,rBSfX) 

CALL  STRLBK(SAR1,IBSAR1,1ESAR1) 

CALL  STRLPM(SAR2,1BSAR2,TKSAR2) 

CALL  STRLEH(SARG,IBSAR3,I£SAR3> 

CALL  STRLEU(SAR4,IBSAR~),':ESAR4) 

CALL  STRLEK(SARS,IBSAR'>,IESAa5) 

CALL  S'n?LER(S«R6,IBSAR6,ILSARC) 

CALL  STRLEN  ( SAR7  ,  IBSAR7 ,  IESAP,7 ) 

C 

isrcn  -  lESARl  +  lESGO 

ISUK2  »  IESAR2  +  7  +IESAR3  +  6 

IS'IHJ  -  ISSAR4  +  IBSH  +  f  +  TZSAR5  +  lESALP 

ISUK4  -  IESAR6  lESA  +  2  +  IESAR7  +  lESFT 

c 

C 

Calculate  the  quantity  ?i»3. 14159267. ,  and  other  paraaeters 
C 

PI-DACOf?(-i,ODOO) 

E.AD3  «  3. •*.5 
G  »  E/(2.*(1>VNEW) ) 

YS  -  Y/RAD3 
OMEGA  “  2.*PI*FREQ 
PRUTT  PI-', PI,'  OHSGAc', OMEGA 

HSTEPS  »  NC'YCLS*KPPCYC 
PRIHT  NSTrFS-',RSTEPE 

PERIOD  -  1,/FREQ 
DELT  -  PERIOD/HPFCyC 
DELGAH  «  (GAMPK2-<AMPK1)/HIHC 
ISTART  ~  HPPCyC/4.  +  1 
lEKD  -  ISTART  +  HFPCYC-l 
C 

Sat  Up  the  loop  for  carrying  out  integration  for  each  step  ia 
C  peak  strain. 

n 

HIHCPl  -  HIKC  +  1 
DO  777?  J-=1,HI1ICP1 
C 

c— —TO-  Clear  the  arrays  for  strets,  strain,  strain  rate,  and  tlxse. 

C 


onnoo  n  nn  uno  ono  one 


nc  10  K*=1,HST2PS 
TAn(K5=0.0 
10  GAKMA(K)«0.0 

“«>«=■>  Update  the  peak  Btrein  and  reinitialize  tijcc  and  1(1). 

T  -  0.0 
Z(l)  -  0.0 

AKP  GAMPK*.  +  EBIiGAM*(J-l) 

PEAKSr{J)  "  AMP 
PEEQST<J)  -  RAD3»AMP.'2. 

•*■“«>  Carry  tat  the  integration  far  the  current  peak  strain 

DO  30  K  -  1, ESTEPS 

CATi  IHTFDH(Z,T,DEI,T,NSQ) 

SAC<K+1)  -  2{1) 

GAMM^(K+1)  «  AMP  •  DSIH(WEGA*T) 

30  CONTUTOE 

«““<►  Cospute  daxiping  for  the  curtoL.t  peak  strain. 

TAUMIH  ~  0.0 
TAUMAZ  «  0.0 

CAM.  HAXIM(TAO,HSTEPS,TADMIH,TAUKAX) 

BNGAB5  »  0. 


DO  50  I-lSTAR‘f,IEHD 
EKGABS  »  EHGABS  + 

> 

50  COHTIHDE 


.5  *  {TAU(1+1)  +  TAU{I))  * 

(GAMKA(1+H  -  GAMMA(I)» 


ENGSTO  -  .5»TAUMAX*AHP 
ENRGAB(J)  SSGABS 
SHRGST{J)  -  EKGSTO 
2TA(J)  »  EKGArS/(a.«PI*EKGSTO) 

WRIT2(«,  ’  (.!X,13,3(3X, 210.4) )  ’)J,P2AJCST(Cf),EPRGAB(J),ETA(J) 
.'777  COHTIHUE 


CAM,  MAXIM(ETA,HIHCPl,EtAMIH,ETA«AX) 

.«->  Hrite  results  to  output  data  file  (unit  10,  12,  13,  and  14)  and  to 
output  text  file  (unit  11). 


WRITE(10,  '  (1X,2(E16.8,  '  '  ' )  j  ‘ )  (PEA31ST(  J)  ,ETA(  J) ,  J-1,HIHCP1 ) 

WRITEdl, '(SX,  ■ '0  '-,12, . ,A,A, . )’) 

>  ISUMl,SARi,SGO/I5SGOjIESGO) 

WRITE!  11,  '  (5X,  •  •!  •  ',12, . ,A,F5,1,2X,A,S6.1, . )  ’  ) 

>  IS'v3H2,SAR2,YS,SAR3,G 

WRITE(11, '(SX,  • ‘2  ’'.12,"  “’‘,A,A,2X,A,A, . )') 

>  ISOM3 ,  EhP.« ,  SN  ( IBSM;  lESH) ,  SAR5,  SALF(  IBSAU?  J  lESALP ) 

WRITE' 11,  •  (5X,  •  *3  'M2, . ,A,A,2X,A,A,  ••■'•)•) 

>  ISUH4 ,  SAR6 ,  SA  ( IBSA  J  lESA ) ,  SAR7 ,  SFT  ( IBSPT  J  lESPT ) 

WRITE{il,  ■  (5X,  "4  38  “gMin-"  ,E9 .4,2Z,  "gHax»"  ,E9 . 4 , . )') 

>  fiTAMIH,ETAHAX 

WRITE!  12,  ■  (1X,2(E16,8,  ","))')  <PBAKST|J)  ,ENRGAB(j;  ,J«<1,KI»CP1) 
WRITE!  13,  ■  (1Z,2!E16.8,  ",’*))•)  (PKEQST{J) , ETA!  J) ,  J»i,HIllCPl ) 
WRl‘tB{14,  •!1X,2!E16.8,  '  ‘ ) )  ' )  (PKEQST(  J)  ,ENRGAB(  J) ,  J-1,HIKCP1) 

PRIHT  MAX.  LOSS  FACTOR- *, ETAHAZ 


1000 

1001 

1200 

1201 

1500 

1501 


FORHAT(A70) 

FORMAT ( 2 AS 0) 

FORMAT!'  ',3X, 'l',5X,I7) 

FORMAT!'  •,3X, '2' .5X,I7) 

FORMAT!'  ■,5X,E10.4,5X,E10.4) 

FORMAT! ■  ’,SX,E10.4,5X,E10.4,5X.E10.4) 

STOP 

S»D 


SUBROUTIHE  DERIV!3,T,ZDOT) 

IMPLICIT  PXAL»8  !A-H,0-Z) 

P.EAL»8  OMEGA,T,Z!l),ZDOT!l),HN 

COMHOK/BLOXl/AMP, OMEGA,  YS,S,G,AI,rA,HH,GAMHA0,A,FT,RAD3 

GAfe»  -  AMP»DSIK!CSffiGA«T) 

GAMBOT  -  AMP*(»iEa\*DCOS(OMEGA»T) 

DLTAGA  o  DABS!GA.MMA)  -  SAM34AC 
IF ! DLTAGA.LT. 0.0) THEN 
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V  o 


B£TA-(l./3.)*(K*MJfA>*(  gamma  -  Z(l)/G  + 

>  Rai)3'’?l'»ERI’(A*GAHHA/RAD3)'‘aMiT(-GAMM&*GaKDOT)  • 
IF{HH.SQ.1.0)THEH 

2DOT(1)  -  G*(  GJaiBOT  -  DAfia(GAMDOT> » (Z (i }-BS?A> /YS  ) 
BX^B 

ZDOt(l)  -  G*(  GAMDOT  -  CABS(GAM)DOT)» 

>  (IJAaS(B(i)-B2'fA>/YS)«»(NH-l)* 

>  {2(l)-BSrA)/YS  ) 

EHDir 

SIfSB 

Zr>M:(l)  »  6«6AM1>OT 
K?5DIF 
REYUHH 
EKD 

>ni«i<w»»K>sDSE .  FOR 


EXAMPLE  OF  AN  m?IJT  Fn.E: 


6HA  Shenr  Bysteretic  Modal  with  Elasticity  Outside  gO 
494  psi 
10.36x106  psi 
3 

.81& 

28.5E06  0.35 

855. 


t 


APPENDIX  F 


FORTRAN  MAINPRQGRAM  (all  subroutines  and  sanction  subprograms  that  are  not 
specifically  related  to  the  calculations  being  made  here  are  given  in  Appendix  H): 


Progran  Naset  BBKDR.rOR 
C 

C*— Thif  program  will  carry  out  o  Runga-Kutta  Istegration  on  the 
C  3tr0»3-t>traiR  aquations  of  the  tiH&  hysteretic  zsodel . 

C  Note  that  inelastic  behavior  is  supressed  beyond  BPSLHO. 

C  "ihm  physical  g«oC;atry  is  ttat  of  a  rectangular  bean  under  bending. 

C  Plane  sections  are  aaattcuo  assuaed  to  reaain  plane. 

C  We  begin  by  specifying  a  B>a::iauii  peak  strain  at  tlio  surface. 

C  A  linear  atrein  distribution  is  assnued  to  exist  about  the 

C  neutral  axis,  tfhe  height  coordinate  (i.a.  y)  is  divided  into  fine 

C  incrensi’ts  and  the  stress  is  cosputed  for  each  of  these  increments. 

C  Based  o,a  the  nonlinear  uniaxial  stress  distribution,  the  torque 

C  is  coaputed  for  each  time  inJtant.  T'he  resulting  curves  for  torque  and 

C  bending  nngle  will  be  stored  in  date  files  and  then  will  be 

C  used  to  calculate  the  loss  factor  (damping)  for  that  pemrticular 

C  peak  strain.  The  peak  strain  is  then  incremented  and  the  process 

C  is  repeated  until  the  final  peak  strain  is  used.  The  results  for 

C  loss  factor  vs.  straiti  are  then  stored  in  an  ASCII  data  file  tc  be 

C  used  in  plotting  with  GRAPHER;  ASCII  data  is  arranged  in  column  pairs. 

C 

c 


c 

c 

c 

c 

c 

c 

c 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

•»» 

A- 

c 

c 

c 

c 

c 

c 

c 

c 


HOKENCLATtlREt 


REAL  COHSTAKTSi 

A  Material  constant  controlling  shape  oZ  hysteresis 

ALFA  “Ey/(E-r.y)  where  Ey  is  t.ie  inelastic  modulus 
AMP  Amplitude  of  strain  input 
B  young ' s  modulus 

EPSLNO  Strain  beyond  which  iuelajtic  growth  is  supressed 

EPSPKl  Minimum  peak  strain 

EPSPK2  Maximum  peak  strain 

FREQ  Frequency  of  cyclic  strain  i.ipui 

FT  Material  constant  corttolling  sice  of  hysteresis 

HEIGHT  Height  (or  thickness)  of  beam 

LSN  Lsegtb  of  beam 

K  Overstress  power  (controls  sharpuesa  of  trc  isitio:.  to  ineX.) 

VHEH  Poisson's  ratio 

WIDTH  Width  (or  depth)  of  Lean 

X  Stress  where  damping  mechanisms  rre  activated 

yy  Height  coordinate  of  bs&m  ctoss-aaction. 

INTEGER  CDKSTAKTSj 

SOCLE  So.  of  cycles  of  oscillation 

HPPCYC  Ho.  of  points  per  cycle  to  be  used  if,  integration 
KINC  Ho.  of  increments  for  tho  range  of  peak  strain 


CHARACTER  STRINGS! 

FILENAME  string  for  filename  assignment  to  a  FORTRAN  -.nit  number 

SALF  string  for  gsatorial  conetant  ALFA 

SA  string  fer  material  constant  A 

SB  string  for  Young’s  moduluc  E 

SBO  string  for  limiting  strain  EPSLHO 

SEP  striug  for  peak  aurfacs  st  .-ain  applied  to  baam 

SET  string  for  material  constant  FT 

SHEI  strina  for  beam  height 

SLBS  string  for  been  length 

SK  string  for  overstress  power  H 

SWXD  string  for  beam  width 

sy  string  for  stress  where  damping  mechanisms  activate!  Y 

TITLE  descriptive  titlo  for  run 


VARIABLE  ARRAYS! 


ANGLE (X) 
ENRGAB(I) 
EKRGST{I) 
EPSL0S(1,K) 

ETh(I) 

PEARST(Z) 

PKEOST(I) 

STPEGSd^JC) 


Angie  of  bending  curvature  (index  K  represents  tisrt) 
Energy  absorbed  in  one  cycle  (I  represents  position) 
Energy  absorbed  in  one  cycle  (I  repretent-  position'. 
Axial  strain  (index  I  rep.  pos  in  cross- -ection, 
index  K  rep.  ti-wo) 

Ix»6E  factor  (I  represents  position) 

Peak  strain  (index  I  represents  position) 

Peak  equivalent  strain  (index  1  represents  position) 
Axial  stress  (Index  I  rep.  pos  in  crocs-seotion, 
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index  K  rep.  tijoe) 

TORQnE(K)  End  torque  acting  on  bean  crosa-section  (K  rep.  tiiie) 

ZjK)  Axial  stress  passed  iron  INTPUH  (index  K  rep.  ti»e) 


SUBRODTIHES: 

DERIV 

INTFHH 

MAXIM 

STRLEK 


Contains  the  differential  eqs.  (invoked  by  IKTPUN) 
Integration  routine  (4th  order  forward  Runge-Kutta) 
Detemines  the  Bin.  £  B^tx.  values  of  an  array 
Counts  the  nuBber  of  characters  in  a  string 


FONCTIOK  SUBPROGRAMS: 

ERF(X)  Error  function  of  X 

FACT(K)  K  factorial 

SGN(V)  SignuB  function  of  V 

UNIT(X)  Unit  step  function  of  X 


IMPLI'’!'^  REAL*8  (A-H.O-X) 

IHTEGER  NEQ 

REAL*8  DERIV,PLOAT,T,TEND,2(20),N,LEH 
DIHEHSION  EPSI>ON(101,1001), 

>  STRESS (101,1001), TORQUE ( 1001 ) , ANGLE ( 1001 ) , 

>  yY(101),ETA(200),PEAKST(200),EKGA(200),ENGS(200) 
CBARACTER*20  FILENAME 

CHARACTER*  70  TZTZiE 

CHARACTER*70  SX ,  SE ,  SH,  SALE ,  SEO ,  SA,  _  ^'T,  SLEN,  SHEI ,  SWID,  SEP 
CHARACTER*  (•)  SARI ,  SARP , SAR2 ,  SAR3 , SAP.4 , SAR5 ,  SAR6 ,  SAR7 ,  SAR8 ,  SAR9 , 

>  SAIO 


parameter  (SARl-‘e0»‘ ) 
PARAMETER  (SARP-'eP-‘) 
PARAMETER  (SAR2-‘y-') 
PARAMETER  ( SAR3~  ‘  E- '  ) 
PARAMETER  '  SAR4- '  n- '  ) 
PARAMETER  (SARS-‘alfa-‘ ) 
PARAMETER  ( SARS-  ‘  a-  ‘  ) 
PARAMETER  ( SAR7- '  f  T- '  ) 
PARAMETER  (SARS-' Length- ' 
PARAMETER  (SAR9-' Height-' 
PARAMETER  ( SAIO-' Width- ' ) 


COMMON/ BLOK 1 /AMP , OMEGA , Y , E , ALFA , N , EPSLNO , A,  FT ,  RATIO 


NEQ  -  1 

— — >  Interactive  input  of  filenemes  for  the  luiterial  data  input  file 
and  tor  output  files. 

WRITE (*,•) 'Enter  the  naae  of  your  input  data  file.' 

IPL-9 

1  READ  (*,2)  FILENAME 

2  FORHAT(A) 

WRITE  (*,  3  )  FILENAME 

3  FORHAT( '  ' ,3X,A20) 

IF(IFL.EQ.9)THEN 

OPEN(IFL, FILE-FILENAME, STATUS-'UNKNOWN'  ) 

ELSE 

OPEN  ( IFL ,  FILE-FILEHAHE ,  STATUS- '  UNKNOWN '  ) 

ENDIF 

IF  (IFL.EQ.9)  THEN 

WRITE (*,*)' Enter  nace  of  file  for  Loss  Fact.  vs.  Peak  strain' 
IFL-10 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.IO)  THEN 

WRITE (*  •) 'Enter  nase  for  tho  summary  file' 

IFL-11 
GO  TO  1 
ENDIF 

IF  (IFL.EQ.ll)  THEN 

WRITE (*,*)' FilenaBe  for  ENGABS  vs.  BNGSTO  at  each  peak  strain' 
IPL-12 
GO  TO  1 
ENDIF 
REWIND  09 
REWIND  10 
REWIND  11 
REWIND  12 
C 

Read  input  quantities  frost  input  file. 

C 

READ(9,5)TITtE 
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RBAD(9, • (A) ' )Sy 
REAI)(9,  •  (A) '  )SE 
READ(9, • (AJ ■ )SH 
READ(9,‘(A)‘)SALP 
R£AD(9,  ■  (A>  •  )SIiN 
REAi>(9,'{A)')SWID 
READ(9,  •(A)')Sa2I 
5  POKMAT(A) 

READ(9,»)E,VHEH 
READ(9,«)i 
READ{9,*)ALI'A 
READ(9,*)H 
KEAD(9,*}LEH 
READ( 9 ,• JWIDTH, HEIGHT 
C 

Interactive  input  of  other  material  parasseterot  strain  eaplitude, 
C  and  nuziber  of  cycles  of  loading  to  be  used  in  calculations. 

C 

WRITE (*,*)•  Enter  EPSIJIO’ 

REAO(5,*)E?SLHO 
WRITE(»,*)*  SPSLHO-'.EPSbHO 

WRITE (*,«)'  Enter  character  string  for  SPSLNO" 

REAO(S, ‘ (A) ‘ )SEO 

WRITS(«, '(XX, ■ ‘SEO-’ 'jA) ‘jSEO 

WRITE(»,»)'  Enter  A  AKD  FT* 

READ(5,*)A,FT 

WRITE(*,»)’  A-',A, ■  FT-', FT 

WRITE (♦,*)'  Enter  character  string  for  a ‘ 

P.EAD(5,  '  (A)  '  )SA 

WPJTE(*,  •  (IX,  •  'SA-'  ’,A)  *  )SA 

WRITE {*,»)■  Enter  character  string  for  fT' 

READ(5,  ■  (A)  •  )SPT 

WRITE(*,  '  (IX,  ‘  'SFiy-'  •  ,A)  •  )SFr 

WRITE (*,*)'  Enter  Ho.  of  points  on  linear  strain  profile* 

READ(5,*)HGRIDP 

WRITE  (*,*)■  HGRIDP- '  ,  KGRIDP 

WRITE (♦,•)■  Enter  FREQ,  EPSPKl  and  BPSPK2  (Kin,  Max  Peak  Strains)’ 
READ ( 5 , * ) FREQ , EPSPKl , EPSPK2 

WRITE(»,«)'  FREQ-', FREQ, •  EPSPKl- ', EPSPKl , '  EPSPK2-’ ,EPSPK2 
WRITE (*,*)'  Enter  character  string  for  EPSPK2 ' 

READ(5, ' (A) ’ )SEP 

write;*, ’ (IX, • 'SEP-'*, A) ’ >S2F 

WRITE (*,♦)'  Enter  Number  of  Cycles  and  Number  of  Points  per  Cycle’ 
READ  (  5 ,  •  )  NCyCLE ,  HPPC'/C 

WRITE (*,»)'  NCYCLE- ' , NCYCLE , ’  HPPCYC- ’ , NPPC YC 

WRITE(*,*)'  Enter  NINC  (No.  of  oncremonts  between  peak  strains)' 
READ(5,*)HINC 
WRITE(*,*)'  NINC-', NINC 
C 

C"— >  Determine  length  of  various  character  strings  for  later  use. 

C 

CALL  STRLEK(SEO,IBSEO,XESEO) 

CALL  STRLEN(SEP,IBSEP,IESEP) 

CALL  STRLEN(SY,IBSY,1ESY) 

CALL  STRLEN(SE,IBSE,1ESE) 

CALL  STRLEN(SN,IBSN,IESH) 

CALL  STRLEN(SALF,IBSALF,IE3ALF) 

CALL  STRLEN(SA,IBSA,I2SA) 

CALL  STRLEN(SFT,IBSFT,IESFT) 

CALL  STRLEN(SLEK,IBSLEN,IESLEH) 

CALL  STRLEN(SHEI,IBSHEI,IESHEI) 

CALL  STRLEN(SWID,IBSWID,IESWIP) 

CALL  STRL£N(SAR1,IBSAR1,IESAR1) 

CALL  strlen;sarf,ibsarp,iesarf) 

CALL  STRLEN(SAR2,IBSAR2,IESAR2) 

CALL  STRLEH(SAR3,IBSAR3,IESAR3) 

CALL  STRLEN(SAR4,IBSAR4,IESAR4) 

CALL  STRLEN(SAR5,IBSAR5,IESAR5) 

CALL  STRLEN(SARe,IBSAK6,XBSAR6) 

CALL  STRLEN(SAR7,IBSAR7,IESAR7) 

CALL  STRLEK(SAR8,IBSAR8,I£SAR8) 

CALL  SjrRLEN(SAR9,IBSAR9,IESAR9) 

CALL  STRLF.N(SA10,IBSA10,IESA10) 

C 

ISUMl  -  lESARl  +  lESEO  +  2  +  lESARP  +  lESEP 

ISUH2  -  IESAR2  +  lESY  +  2  +IESAR3  +  lESE 

ISUM3  »  IESAR4  +  lESN  +  2  +  IESAR5  +  lESALF 

ISt7M4  -  IESAR6  +  lESA  +  2  +  IESAR7  +  lESFT 

ISOMS  -  IESAR8+IESLEN+2+IESAR9+IESHEI+2+IESA10+IESWID 
C 

C— — >  Calculate  the  quantity  Pi-3.14159267,..,  and  other  parameters 
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09*^  non  noon  oooono  non  ooonn 


PI-DACX)S{-1.0D00) 

OMEGA  -  2.*P1*PREQ 

PRINT  P1-',PI,  ■  OMEGA-*, OMEGA 

KSTEPS  -  NCYCXE*NPPCYC 

NSTPPl  -  HSTEPS+i 

PRINT  NSTEPS-’jNSTEPS 

PERIOD  •*  l./FREQ 

DELY  «•  KEIGHT/NGRIDP 

DELT  -  PERIOD/NPPCYC 

DELEPS  -  {EPSPXY-EPfPKlj/NINC 

ISTART  -  NPPCYC/4.  +  1 

lEND  -  ISTART  +  KPPCYC-1 


“->  Set  ap  radial  grid. 

DO  7,  I-X,HGRIDP 

7  YY(I)  -  -{HEIGKT/2.)  +  DELY*(2.*I  -  I.)/2. 

Set  up  the  loop  for  carrying  out  Integration  for  each  iitf-p  in  peedc 
■train  and  position  in  the  cross-section. 


DO  7777  J-1,NIHC+1 

■— >  Reinitialize  T  and  Z(l).  Clear  the  arrays  used  in  the  integration  loop. 

T  -  0. 

Z(l)  -  0. 

DO  9  K-1, NSTPPl 
DO  8  I-1,MGRIDP 
STRESS(I,K)  -  0. 

8  EPSL0N{I,K)  -  0. 

9  CONTINUE 

■“»>  Update  the  peak  surface  strain. 


AMP  -  EPSPK1+DELEPS*(J-1) 

PEAKST(J)  -  AMP 

—« «>  Carry  out  the  integration  at  each  point  in  one  half  of  the  beam 

cross-section  for  the  current  surface  strain  (only  one-half  of  the 
cross-section  need  be  considered  because  of  symmetry  ed>out  the 
neutral  axis ) . 

DO  11  I-l,SGRIDP/2 
T  -  0. 

Z(l)  -  0. 

EPSL0N(I,1)  -  0. 

STRESS(I,1)  -  0. 

RATIO  -  yY( I) /(HEIGHT/2.) 

DO  10  K  -  l.HSTEPS 

CALL  INTFUN(S,T,DELT,NEQ) 

EPSL0N(I,K+1)  -  -AKP  •  DSIN(OHEGA»T)  •  RATIO 
STRESS(I,K+1)  -  2(1) 

10  CONTINUE 

11  CONTINUE 

Using  symmetry,  determine  stresses  and  strains  for  the  other 
half  of  the  cross-section 


DO  13  I-KGRIDP/2+l,NGRIDP 
DO  12  K-1, NSTPPl 

EPSLON(I,K)  -  -EPSLON(NGRIDP-I+l,K) 

12  STRESS(I,K)  -  -STRESS(NGRIDP-I+1,K) 

13  CONTINUE 


“■— >  Compute  the  acting  end  torque  and  angle  of  beam  curvature. 

DO  15  K-1, NSTPPl 
15  TORQUE(K)»0.0 

DO  25  K-1, NSTPPl 
DO  20  I-1,NGRIDP 

20  TORQUE (K) -TORQUE (K)  +  (-YY(I) )«STRESS(I,K)*BELY 

TORQUE (K)  -  WIDTH*TORQUE(K) 

ANGLE (K)  -  2. *LEN*EPSLOH(l,K) /HEIGHT 
25  CONTINUE 

— “>  Compute  damping  for  the  current  pe2dc  surface  strain. 

CALL  MAXIM (TORQUE, NSTPPl, T0RM1N,T0RHAX) 
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CALL  MAXIM(AN6LE,HSTPPl,AN(aiIN,ANGHAX) 

ENGABS  -  0. 

C 

DO  50  I-ISTART.IEHD 

EKGABS  -  ENGABS  +  .5  »  ( TORQOE { I+l )  +  TORQUE(I) )  » 

>  (AKGLR{T-M)  -  ANGLE(I)) 

50  CONTINUE 
C 


ENGSTO  -  .  5*TORMAX*ANGMAX 

ENGA(J)  -  ENGABS 

ENGS(J)  -  ENGSTO 

ETA(J)  -  ENGABS/ ( 2. *PI*ENGSTO) 

WRITE(*, • (3X,I3,2(3X,E10.4)) •)J,PEAKST(J>,ETA(J) 

7777  CONTINUE 

CALL  MAXIH(£TA,HINC,ETAHIK,ETAHAX) 

C 

C«»>«— >  Write  reiult*  to  output  data  file  (unit  10)  and  to 
C  output  text  file  (unit  11). 

c 

HRITE(10,  •(1X,2(E16.8,  '  * ,  ‘  *  )  )  *  )  (PEAKST(  J)  ,ETA(  J)  ,  J-1  ,NIHC) 
WRITE(12, • (1X,2(E16.8, ‘ •, ' ' ) ) ‘ ) (ENGS( J) ,ENGA( J) , J-1,NINC) 
WRITE(11,  ■  (5X,  ■  *  M2, . ,A,A,2X,A,A, . )  ') 

>  ISUHl , SARI , SEO ( IBSEO : lESEO ) , SARP , SEP { IBSEP i lESBP } 

WRITS(11,  ■  (5X,  •  •M2,‘'  '*  •,A,A,2X,A.A, . )‘) 

>  IS0K2,SAR2,Sy(IBSYsIESY),SAR3,SE(IBSEjISSE) 

WRITE(n,  ■  (5X,  '  ‘2  *  •  ,12, . ,A,A,2X,A.A, . )  *  ) 

>  ISUM3 ,  SAR4 ,  SN  ( IBSN  :  lESN  )  ,  SAR5  ,  SAL?  ( IBSALF :  lESALF ) 

WRITE(11,  *  (5X,  •  '3  •  •  ,12, . ,A,A,2X,A,A, . )  ’  ) 

>  iaUM4 ,  SAR6 ,  SA(  IBSA  S  lESA  )  ,  SAR7  ,  SPT  ( IBSPT  :  lESPT ) 

WRITE(11,  •  (5X,  "4  ",I2, . ,A,A, 2X,A,A, 2X, A,A, . )  ‘  ) 

>  ISnH5,SAR8,SLEN(IBSLEH:IESLEH), 

>  SAR9,SHEI(IBSHEl!lESHEI), 

>  SA10,SWID(IBSWIDsIESWID) 

WR1TE{11, '(SX, • '5  26  "Peak  Loflfi  Factor-’ ' ,E9 , 3, . ) ’ )ETAMAX 

PRINT  •,*  MAX  LOSS  FACTOR-’ ,ETAMAX 
C 


1000  FORHAT(A70) 

1001  FORMAT(2A60) 

1200FORMAT(’  ‘ , 3X, ’ 1 ' , 5X, 17 ) 

1201PORMAT(’  ’,3X, '2' ,5X,I7) 

1500FORMAT(’  ' , 5X,E10 .4 , 5X,E10.4 ) 

1501  PORMATC  ’,5X,E10.4,5X,E10.4,5X,E10.4 
STOP 
END 


) 


SUBROUTINE  DERIV(2,T,ZDOT) 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  NEQ 

REAL*8  C»ffiGA,T,Z(20),ZDOT(20),NN 

COHHON/BLOKl  /AMP ,  OMEGA,  Y ,  E ,  ALFA,  HH  ,  EPSLNO  ,  A .  FT,  RATIO 
C 

EPSLON  »  -AMP  *  DSIN(OMEGA*T)  •  RATIO 
EPDOT  -  -AMP  ♦  OMEGA  •  DCOS(OHEGA»T)  •  RATIO 
DLTAEP  -  DABS (EPSLON)  -  EPSLNO 
IF ( DLTAEP . LT . 0 . 0 ) THEN 

BETA-(E*ALFA)*(  EPSLON  -  Z(l)/E  + 

>  FT*ERF(A*EPSLON)*UNIT(-EPSLON*EPDOr)  ) 
IP(HH.EQ.1.0)THEH 

ZDOT(l)  -  E*(  EPDOT  - 

>  DABS(EPDOT)»(Z(l)-BETA)/Y  ) 

ELSE 

ZDOT(l)  -  E®(  EPDOT  -  DABS(EPDOT)* 

>  (DABS(Z(l)-BETA)/Y)«»(NN-?.)» 

>  (2(1)-BETA)/Y  ) 

ENDIF 

ELSE 

ZDOT(l)  -  B'EPDOT 
ENDIF 
RETURN 
END 
C 

C— — >BENDR.  FOR 


EXAMPLE  OF  AN  INPUT  HLE: 


SMA  Hysteretic  Model  with  Elasticity  Outside  eO 
855  psi 
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28.5x106  psi 
3 

.818 
20  in 
1  in 
.5  in 

28.5E0e  0.35 

855. 

0.818 

3 

20. 

1.  .5 


EXAMPLE  OF  A  BATCH  FILE  TO  RUN  BENDR  ON  THE  CODE  281  MICROVAX: 


set  def  [.hysteresis] 

r  bendr 

bendr.inp 

bendr.dat 

bendr.txt 

benrg.dat 

l.e-4 

l.Oe-4 

80000.  4.e-5 

80000 
4.0e-S 
20 

1.0  .25e-4  1.753-4 

1.75e-4 
2  100 
20 
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APPENDIX  G 


FORTRAN  MAIN  PROGRAM  (all  subroutines  and  function  subprograms  that  are  not 
specifically  related  to  the  calculations  being  made  here  are  given  in  Appendix  H): 


Program  Name)  TORRANGE.POR 
C 

This  program  will  carry  out  a  Runge-Xutta  integration  on  the 
C  stress-strain  equations  of  the  SHX  hysteretic  model. 

C  Hote  that  inelastic  behavior  is  supressed  beyond  GhHKAO. 

C  The  phyuical  geometry  is  that  of  a  solid  cylinder  under  torsion. 

C  Plane  stress  conditions  are  assumed  (pure  sheu  to  be  specific). 

C  We  begin  by  specifying  a  minimum  peak  strain  at  the  surface. 

C  A  lineeu:  strain  distribution  is  assumed  to  exist  in  the  circular 

C  cross  section.  The  radial  coordinate  is  divided  into  fine 

C  increxients  2md  the  stress  is  computed  for  each  of  these  increments. 

C  Zased  on  the  shear  stress  distribution,  the  resulting  torque  is 

C  CKiputed  for  each  time  instant.  The  resulting  curves  for  torque  and 

C  overall  sample  twist  will  be  stored  in  data  files  and  then  will  be 

C  used  to  calculate  the  loss  factor  (damping)  for  that  particular 

C  peak  strain.  The  peak  strain  is  then  increiDented  and  the  process 

C  is  repeated  until  the  final  peak  strain  is  used.  The  results  for 

C  loss  factor  vs.  strain  are  then  stored  in  an  output  file  to  be  used 

C  in  platting  with  GRT^HER;  ASCII  data  is  arranged  in  column  pairs. 

C 

C 

C  NOHENCLATUREs 
C 

C  REAL  CONSTANTS: 

C  A  Material  constant  controlling  shape  of  hysteresis 

C  ALFA  »Ey/(E-Ey)  where  Ey  is  the  axial  inelastic  modulus 

C  AMP  Amplitude  of  strain  input 

C  OIA  Diameter  of  shaft 

C  E  Young's  modulus 

C  G  Shear  modulus  {-E/[2( 1+VNEW) ]} 

C  GAMMAO  Strain  beyond  which  inelastic  growth  is  supressed 

C  GAMPKl  Minimum  peak  strain 

C  GAKPK2  Maximum  peak  strain 

C  FREQ  Frequency  of  cyclic  strain  input 

C  FT  Material  constant  controlling  size  of  hysteresis 

C  LEN  Length  of  Shaft 

C  K  Overstress  power  (controls  sharpness  of  transition  to  inel.) 

C  RADIUS  Radius  of  shaft 

C  VNEW  Poisson's  ratio 

C  Y  Axial  stress  where  damping  aeciiani&ms  ore  activated 

C  YS  Shear  stress  wliere  oamping  mechanisms  are  activated 

C 

C  INTEGER  CONSTANTS  I 

C  NCYCLE  Ho.  of  cycles  of  oscillation 

C  NFPCYC  Ho.  of  points  per  cycle  to  be  need  in  integration 

C  NINC  No.  of  increments  for  the  range  of  peak  strain 

C 

C  CIWRACTER  STRINGS: 

C  PaLSKAME  string  for  fileraraa  assigiusont  tc  a  FORTRAN  unit  number 

C  SALF  string  for  material  constant  *J:.FA 

c  SA  string  for  .aaterial  constant  A 

C  SDIA  string  for  geometric  constant  DIA 

C  SG  string  for  Yctng's  mociuius  G 

C  SGO  string  for  limiting  strain  GAMMAO 

C  SPT  string  i'or  material  con8te:;t  FT 

C  SLEN  string  for  geometric  constant  LEN 

C  SH  string  to^-  overstress  pover  S 

C  SYS  iftring  for  shear  stress  where  damping  mechanisms  activate 

C  TITLE  descriptive  tit’.e  for  run 

C 

C  VARIABLE  ARRAYS; 

C  ANGLE  (A)  Angle  of  twist  of  shaft  (i.-ids:.  Z  represents  time) 

C  £BRGA3(I)  Energy  absorbed  in  one  cycle  (I  repreaer.ts  position) 

C  SBRGST(I)  Energy  absorbed  is  one  cycle  (I  represents  positioni 

C  GAjet-  I,R)  Shear  strain  (index  i  rep.  position,  R  rep.  time) 

C  ETA  Loss  factor  (1  represents  position) 

C  PEAKsr{I)  Peak  strain  (index  I  represent g  position) 

C  PKEQST(X!  Peak  equivalent  strain  (index  I  rspraserts  position) 

C  TAU(I,K)  shear  stresj  (index  I  rap.  positi.oi.,  S  rep.  time) 

C  TORQUE (K)  Resultant  tcrque  acting  on  sliaft 

C  Z(K)  Shear  ctrass  passed  fres'  IKTIUN  (index  K  rep.  time) 
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oonn  n  o n n n o o n n n o n n n n n 


SUBROUTINES i 
DBRIV 
INTPUH 
KAXIH 
STRLEN 


Contains  the  differential  eqs.  (invoked  by  IHTFUN) 
Integration  routine  (4th  order  forward  Runge-Kutta) 
Determines  the  min.  t  moc.  values  of  an  array 
Counts  the  number  of  characters  in  a  string 


FUNCTION  SUBPROGRJVHS: 

ERF(X/  Error  function  of  X 

FACT(K)  K  factorial 

SGN(V)  Signum  function  of  V 

UNIT(X)  Unit  step  function  of  X 


IMPLICIT  REAL* 8  (A-B,0-Z) 

INTEGER  NEQ 

REAL*8  DERIV, FLOAT, T, TEND, TOL,Z(20),N,LEN 
DIMENSION  GAMMA(101, 1001) ,TAU(101, 1001), 

>  TORQUE ( 1001 ),TWIST( 1001), 

>  R(lOl)  ,ETA(200)  ,PEAKST{200)  ,ENGA(200)  ,EHGS(200) 
CBARACTE?*20  FILENAME 

CHARACTER*  7  0  TITLE , S YS , SG , SN , SALF , SGO , SA, SFT , SLEN , SDIA, SGP 
CHARACTER* (*)  SARI , SAR2 , SAR3 , SAR4 , SAR5 , SAR6 , SAR7 , SAR8 , SAR9 , SARP 
PARAMETER  ( SARI-  ‘  gO- '  ) 

PARAMETER  ( SARP-  ’  gP- '  ) 

PARAMETER  (SAR2-'Y8-*) 

PARAMETER  (SAR3-‘G-') 

PARAMETER  (SAR4-'n-*) 

PARAMETER  ( SARS- '  a- '  ) 

PARAMETER  ( SAR6- '  a- '  ) 

PARAMETER  ( SAR7-  ‘  f  T- '  ) 

PARAMETER  ( SARS- '  Length- '  ) 

PARAMETER  ( SAR9- ' Diameter- * ) 

COHHON/BLOKl/AMP ,  OMEGA,  YS ,  E ,  G,  ALFA,  N,  GAMMAO,  A,  FT,  RAD3 ,  RATIO 
K2Q  -  1 

interactve  input  of  filenames  for  the  material  data  input  file 
and  for  output  files. 

WRITE (*,*)' Enter  the  name  of  your  input  data  file.' 

IFL-9 

1  READ (*,2) FILENAME 

2  FORMAT (A) 

WRITE (*, 3 ) FILENAME 
'l  FCSMATC  •,3X,A20) 

IF(IFL.EQ.9)THEN 

OPES(IFL,FILE-FILENAME,STATUS-'UNKNOKN'  ) 

ELSE 

OPEN ( IFL, FILE-FILENAME , STATUS- ' UNKNOWN ' ) 

ENDIF 

IF  (IFL,EQ.9)  THEN 

WRITE(',*) 'En^er  name  of  file  for  Loss  Fact.  vs.  Peak  Strain’ 
IF.’,-10 
GO  Tu  1 
ENDIF 

IF  (IFL.EQ.IO)  THEN 

WRIIE)*,*) 'Enter  name  for  the  summary  file’ 

j;fl*-ii 

GO  TO  1 
ENDIF 

IF  (IFI.EO.ll)  THEN 

WRITE(*,*) 'Filename  for  ENGABS  and  ERGSTO  at  each  peak  strain.' 
IPL-12 
GO  TO  1 
ENDI? 

REWIND  09 
REWIND  10 
SEWTNO  21 
REUIKD  12 
C 

Read  input  quantities  fr<so  input  file. 

READ (9, 5) TITLE 
READ(9  ' (A) ’ )SYS 
REiJ)(9,  '  (A)  •  )SG 
READ(®,  '(A)  ’)SN 
READ(9, ■ (A) •)SALP 
a'JAD(9,  '  (A)  ’  )SLEH 
REA&{9|'{a}‘>SDIA 
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5  FORHftT(&) 

READ{9,*)E,VKEW 
READ(9,*)Y 
READ( 9 , * )ALFA 
READ(9,*)N 
REA0(9,*)LEN 
REAO(9,*)DIA 
C 

C“-~>  Interactive  input  of  other  naterial  paraiiAters,  strain  amplitude, 
C  and  nunber  of  cycles  of  loading  to  be  used  in  calculations. 

C 

WRITE(»,*)'  Enter  GAMMAO* 

REAO(5,*)GAKHAO 
WRITB(*,*)'  GAMMAO-*  .GAMMAO 

WRITE (*,*)'  Enter  character  string  for  GAMMAO ' 

READ(5, ' (A) ' )SGO 

WRITE(*, • (IX, * 'SGO-' • ,A) ’ >SG0 

WRITE(*,*)'  Enter  A  AND  FT’ 

READ (5,*) A, FT 

WRITE(*,*)*  A»',A, *  FT-’, FT 

WRITE (*,•)’  Enter  character  string  for  a’ 

READ(5, ’ (A) ’ )SA 

WRITE(*,  ’  (IX,  ’  •SA-”,A)  ’  )SA 

WRITE (*,*)’  Enter  character  string  for  fT’ 

READ(5, ’ (A) ’ )SPT 

WR1TE(«, ’ (IX, ’ ’SFT-’ ’.A) ’ )SFT 

WRITE(»,*)‘  Enter  No.  of  points  on  linear  strain  profile’ 

READ(5,*)NGRIDP 

WRITE (•,*)'  NGRIDP- ’ , NGRIDP 

WRITE(»,*)’  Enter  FREQ,  GAMPKl  and  GAMPK2  (Min,  Max  Peak  Strains)’ 
READ  ( 5 ,  *  )  FREQ  ,  GAMPKl ,  GAMPE2 

WR1TE(»,*)’  PREQ-’ .FREQ,  ■  GAMPKl-’ , GAMPKl,  ’  GAMPK2-’ ,GAMPK2 
WRITE (*,*)'  Enter  character  string  for  GAMHAP ’ 

R£AD(5, ’ (A) ‘ )SGP 

WRITE(*, ’ (IX, ' 'SGP-’ ’ ,A) ■ )SGP 

write!*,*) ‘  Enter  Nunber  of  Cycles  and  Nuaber  of  Points  per  Cycle* 
READ ( 5 , * ) NC YCLE , NPPC YC 

WRITE (*,*)’  NCYCLE- * , NCYCLE , ’  NPPCYC- ’ , NPPCYC 

WRITE(*,*)'  Enter  NINC  (Ho.  of  increnents  between  peak  strains)’ 
READ(5,*)NINC 
WRITE(*,»)*  NINC-', NINC 
GAMHAP  -  GAMPK2 
C 

Cran>:>  Determine  length  of  various  character  strings  for  later  use. 

C 

CALL  STRLEN(SGO,IBSGO,XESGO) 

CALL  STRLEN(SGP,IBSGP,IESGF) 

CALL  STRLEN(SYS,IBSYS,IESYS) 

CALL  STRLEN(SG,IBSG,IESS) 

CALL  STRLEN(SH,IBSN,IESH) 

CALL  STRLEN(SALF,IBSALF,IESALF) 

CALL  STRLEK(SA,IBSA,IESA) 

CALL  STRLEN(SFT,IBSFT,IESFT) 

CALL  STRLEN(SLEN,IBSLEH,IESLEN) 

CALL  STRLEN(SDIA,IBSDXA,IESDIA) 

CALL  STRLSN(SAR1,IBSAR1,IESAR1) 

CALL  STKLEN(SARP,IBSARP,IESARF) 

CALL  STRLEH(SAR2,IBSAR2,IESAR2) 

CALL  STRLEN(SAR3,IBSAR3,IESAR3) 

CALL  STRLEN(SAR4,IBSAR4,IESAR4) 

CALL  STRLEN(SAR5,IBSAR5,IESARS) 

CALL  STRLEN(SAR6,IBSAR6,IESAR6) 

CALL  STRLEN(SAR7,IBSAR7,IESAR7) 

CALL  STRLEN(SAR8,IBSAR8,IESAR8) 

CALL  STRLEH(SAR9,IBSAR9,I£SAR9) 

C 

ISOMl  -  lESARl  +  lESGO  +2  +  lESARP  +  lESGP 

ISUM2  -  IESAR2  +  lESYS  2  +  IESAR3  +  lESG 

IS0M3  -  IESAR4  +  lESN  +  2  +  IESAR5  +  lESALF 

ISUM4  -  IESAR6  +  ISSA  +  2  +  IESAR7  +  lESFT 

ISUM5  -  I2SAR8  +  lESLEN  +  2  +  IESAR9  +  IBSDIA 

C 

C— — >  Calculate  the  quantity  Pi-3.14159267...,  and  other  parameters 
C 

PI  -  DACOS(-l.ODOO) 

RAD3  -  3. **.5 
G  -  E/(2.*(1+VNEW) ) 

YS  -  Y/RAD3 

OMEGA  -  2.*PI*FREQ 

PRINT  *,’  PI-’, PI,’  OMEGA-’, OMEGA 

NSTEPS  -  NCYCLE*NPPCYC 


DTRC-SME-91/34 


59 


HSTPPi  -  NSTEPS+1 
PRINT  HSTEPS-',HSTEPS 
PERIOD  -  l./FREQ 
RADIUS  >  DIA/2. 

DELR  -  (DIA/2.  )/NGRIDP 
DELT  -  PERIOD/NPPCYC 
PRINT  •  ,  ■  DELR-  *  ,  DELR,  '  DELT-  ’  ,DELT 
DELGAM  -  (GAHPK2-GAMPK1)/NINC 
PRINT  *,  'DELGAM-*  ,DEICAM 
ISTART  -  NPPCyC/4.  +  1 
lEND  -  ISTART  +  NPPCYC-l 
C 

C— — >  Set  up  radial  grid, 

c 

DO  7,  I-1,HGRIDP-1 

7  R(I)-DELR*1 

R  ( NGRIDP ) -RADIUS 
C 

C— — >  Set  up  the  loop  for  carrying  out  integration  for  each  step  in 
C  peak  strain. 

C 

DO  7777  J-1,NINC+1 
C 

C— — >  Cleetr  the  arrays  for  shear  stress  and  strain. 

C 

DO  9  K-X,NSTPP1 
DO  8  I-l, NGRIDP 
TAU(I,K)  -  0. 

8  GAMMA(I,K)  -  0. 

9  CONTINUE 
C 

C— — >  Update  the  peak  strain  and  reinitialize  tine  and  2(1). 

C 

T  -  0. 

2(1)  -  0. 

AMP  -  GAMPKl-f-DELGAM*(J-l) 

PEAKST(J)  -  AMP 


C 

C— — >  Carry  out  the  integration  at  each  point  along  the  radial 
C  coordinate  of  the  circular  cross-section,  for  the  current 

C  surface  strain. 

C 

DO  11  I-l, NGRIDP 
T  -  0. 

2(1)  -  0. 


GAHHA(X,1)  •  0. 

TAU(I,1)  -  0. 

RATIO  -  R( I) /RADIUS 
DO  10  K  -  1,NSTEPS 

CALL  INTPUN(2,T,DELT,NEQ) 

TAU(I,K+1)  -  2(1) 

GAMMA(I,K+1)  -  AMP  •  DSIN(OMEGA»T)  *  RATIO 

10  CONTINUE 

11  CONTINUE 
C 

C-—- >  Conpute  the  acting  torque  emd  angle  of  twist  of  the  shaft, 

c 

DO  15  K-1,NSTPP1 
15  TORQUE(K)-0.0 
C 

DO  25  K-1,HSTPP1 

DO  20  I-l,  HGRIDP-1 

20  TORQUE (K) -TORQUE (K)  +  R(I)*TAU(I,K)»2.*PI*R(I)»DELR 

TORQUE (K) -TORQUE (K)  +  (R( NGRIDP )-DELR/4. )*TAU(NGRIDP,K)* 

>  2.*PI»(R(NGRIDP)-DELR/4. )»DELR/2. 
TWIST(K)  -  GAMMA ( NGRIDP, K)»LEN/RADIUS 

25  CONTINUE 
C 

C— — >  Conpute  danping  for  the  current  peak  strain. 

C 

C 

CALL  MAXIM ( TORQUE, NSTPPl,TORMIN,TORMAX) 

CALL  MAXIM ( TWIST, NSTPF1,THIMIN,TWIMAX) 

ENGABS  -  0. 

C 

DO  50  I-ISTART,IEND 

ENGABS  -  ENGABS  +  .5  •  (TORQDE(I+l)  +  TORQDE(I))  * 

>  (TWIST(I+1)  -  TWIST(I)) 

50  CONTINUE 

C 

ENGSTO  -  .5*TORMAX*TWIMAX 
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EIIGA(J)  >  ENGABS 
ENGS(J)  -  ENGSTO 
STR(J)  -  EHGABS/(2.»PI*ENGSTO) 

WRITB(», • <3X,I3,2(3X,E10.4) ) ‘ )J,PEAKST(J) .ETA(J) 
7777  CONTINUE 

CALL  KAXIH(ETA,NINC,ETAMIN,ETAMAX) 


»>  Write  results  to  output  data  file  (unit  10,  12,  13,  and  14)  and  to 
output  text  file  (unit  11). 


WRITE(10, •(1X,2(E16.8, ' * , ’ ' ) ) ‘ ) (PEAKST( J) ,ETA( J) , J-1,NINC) 
WRITE(12, •(1X,2(E16.8, •■,•*))') (2MGS(J) ,ENGA{J) ,J“1,NINC) 
WRITE(11,  •  (5X,  •  ’0  ”,I2, . ,A,A,2X,A,A, . )  '  ) 

>  ISUKl , SARI , SGO ( IBSGO : lESGO ) , SARP , SGP ( IBSGP ! TESGP ) 

WRITE ( 11,  '(SX,  "I  ",12, . ,A,A,2X,A,A, . )  ‘  ) 

>  ISDM2 , SAR2 , SYS ( IBS YS  t lESYS ) , SAR3 , SG ( IBSG : lESG ) 

WRITE(11,  ■  (5X,  "2  ■  ',12, . ,A,A,2X,A,A, . )  ‘ ) 

>  ISUH3,SAR4,SN(IBSN:IBSN),SAR3,SALE(IBSALP:IESALF) 

WRITE(11,  ■  (5X,  "3  "  ,12, . ,A,A,2X,A,A, . )  ‘  ) 

>  ism44 , SAR6 , SA( IBSAt lESA) ,SAR7 , SPT( IBSPT: lESPT) 

WRITE ( 11,  ‘(BX,  "4  ",I2, . ,A,A,2X,A,A, . )  *  ) 

>  ISUM5 .  SAR8 ,  SLEN  ( IBSLEN  t  lESLEN  )  ,  SAR9  ,  SDIA  ( IBSDIA :  lESDIA ) 

WRITE(11, ’(SX,  "5  26  'g- "  ,E9 .4 ,2X,  ‘ ’SDC-"  ,E9 . 4 , . )‘)ETA,SDC 

WRITE(11, ‘(5X, ■ '6  32  ‘Peak  Loss  Factor-’ ‘ ,E9. 3, . ) ’ )ETAMAX 

PRINT  MAX  LOSS  FACTOR- ‘ , ETAMAX 


1000 

1001 

1200 

1201 

1500 

1501 


FORMAT  (A70) 

FORMAT  ( 2  A60) 

FORMATC  •  ,3X, ’1’ ,5X,I7) 

FORMAT ( ■  * ,3X, *2 ■ ,5X,I7) 

FORMATC  ',5X,E10.4,5X,E1Q.4) 

FORMAT( ■  •,5X,E10.4,5X,E10.4,5X,E10.4) 

STOP 

END 


SUBROUTINE  DERIV(Z,T,ZDOT) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

INTEGER  NEQ 

REAL*8  OMEGA, T,Z(20),ZDOT(20),NH 

COMMON/BLOK  1 /AMP ,  OMEGA ,  YS ,  E  ,  G ,  ALFA ,  NN  ,  GAMHAO ,  A ,  FT ,  RAD3 ,  RATIO 

GAMMA  -  AMP  *  DSIN(OHEGA*T)  *  RATIO 
GAMDOT  -  AMP  »  OMEGA  •  DCOS(OMEGA*T)  •  RATIO 
DLTAGA  -  DABS  (GAMMA)  -  GAMHAO 
IF ( DLTAGA . LT . 0 . 0 ) THEN 

BETA-(l./3.  )*(E*ALFA)*(  GAMMA  -  Z(l)/G  + 

>  RAD3*FT*ERF(A*GAMMA/RAD3)*UHIT(-GAMMA*GAMDOT)  ) 
IF(NN.EQ.1.0)THEN 

ZDOT(l)  -  G»(  '^AMDOT  -  DABS(GAMDOT) »(Z ( 1  )-BETA) /YS  ) 
ELSE 

ZDOT(l)  -  G*(  GAMDOT  -  DABS(GAHDOT)* 

>  (DABS(Z(1)-BETA)/YS)**(NN-1)* 

>  (Z(1)-BETA)/YS  ) 

ENDIF 

ELSE 

ZDOT(l)  -  G*GAMDOT 
ENDIF 
RETURN 
END 


— >TORRANGE . FOR 

EXAMPLE  OF  AN  INPUT  RLE: 


SHA  Hysteretic  Shear  Model  with  Elasticity  Outside  gO 
494  psi 
10.56x106  psi 
3 

.818 
2  in 
.5  in 

28.5E06  0.35 

855. 

0.818 

3 

2. 

.5 
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EXAMPLE  OF  A  BATCH  FILE  TO  RUN  TGRRANGE  ON  THE  CODS  281 MICRQVAX: 


set  def  [ .bysteresis] 

r  torrange 

torrar.ge .  inp 

torrange.dat 

torrange.txt 

torrenrg.dat 

1.56e-4 

1.56e-4 

eOOOO.  4.a-5 

80000 

4.0e-5 

40 

1.0  .05e-4  2.57e-4 

2.57e-4 
2  100 
20 
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APPENDIX  H 


GENERAL  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  USED  IN  THE  FORTRAN 
PROGRAMS  IN  APPENDICIES  B-G: 


c  title:  strlen.for 


c 

c  author:  shw 

c  Turbulence  Research  Laboratory 

c  SDNY  «  Buffalo 

c 

c  trlscottgubvKs 

c 

c  function: 

c  find  the  first  and  last  cheiracter  in  a  string 

c 

c  inputs: 

c  string  -  character  string 

c 

c  outputs: 

c  ib  -  first  non-blank  character 

c  ie  -  last  non-blank  character 

c 

c  subroutines  required: 
c 

c  nodif icaticns : 

c 

c 

subroutine  strlen( string,  ib,  ie) 
character  string* (*) 
j  ••  len( string) 
do  10  i=j,l,-l 

if (string(i:l)  .ne.  ‘  ')  goto  20 
10  continue 

ie  ■«  1 
ib  -  1 
return 
20  ie  i 

do  30  i“l,j 

if (string(l:l)  .ne.  ‘  ')  goto  40 
30  continue 

40  ib  -  i 

return 
end 


SUBROUTINE  MAXIM {r,HPOINT,PHIK,FHAX) 
IMPLICIT  REAL* 8  (A-H,0-Z) 

UIMENSION  F(lOOl) 

FMAX  -  0. 

FKIK  -  0. 

DO  10  I-l,NPOIHT-l 

IF{F{I+1).GT.FMAX)THEN 
FMAX  -  F(I+1) 

ENDIF 

IF ( F ( I+l ) . LT. FMIN ) THEN 
FMIH  -  F(I+1) 

ENDIF 

10  CONTINUE 
RETURN 
END 


REAL* 8  FUNCTION  ONIT(X) 
REAL*8  X 
IF(X.LT.O.0)THEN 
UNIT-0.0 
ELSE 

UNIT-1.0 

ENDIF 

RETURN 

END 
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o  o  o  o 


REAL*8  FDHCTIOH  ERF(X) 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  K 

PI  -  DACOS(-l.ODOO) 

IF(X.GT.3.0D00)THEH 
ERP  -  l.ODOO 
GO  TO  500 

ELSEIF ( X . LT . -3 . ODOO ) THEN 
ERP  -l.ODOO 
GO  TO  500 
ELSE 

TOL  -  .0001 
K  »  0 
KMAX  -  35 
SERIES  •>  1. 

THOLT  -  2.*X/(PI**.5) 

100  K  »  K  +  2 

IF(X.GT.EH10C)THEt! 

PRINT  K  -  KKAX  ERP(  *  ,X»  ' 

GO  TO  500 
ENDIP 

PRIOR  -  SERIES 

ODDTRM  -  (-!)*(  DABS(X)**(2.*{K-1) )  )/ 

>  (  PACT{X-1)*(1+2»(K-1) )  ) 

EVHTRM  -  (  DABS(X)**(2.*K)  )/(  FACT(K)*(1+2*K)  ) 
SERIES  -  SERIES  -i-  ODDTRM  +  EVHTRM 
ERFl  -  THULT  *  PRIOR 
ERF2  “  TMULT  *  SERIES 
DELTA  -  ERP2  -  ERFl 
ADBLTA  -  ABS( DELTA) 

IF ( ADELTA . LB . TOL ) THEN 
ERF  -  ERF2 
GO  TO  500 
ELSE 

GO  TO  100 
ENDIF 
ENDIP 
500  RETURN 
END 


REAL* 8  FUNCTION  FACT(K) 
IMPLICIT  REAL*8  (A-H,0-Z) 
INTEGER  I,K 

IF(K.EQ.O.OR.K.EQ.1)THEN 
FACT  -  1.0 
GO  TO  20 
ELSE 

FACT  “  1.0 
DO  10, 1-2, K 
FACT  -  FACT* I 
10  CONTINUE 
ENDIF 
20  RETURN 
END 


RBAL»8  FUNCTION  SGN(V) 
REAL*8  V 
IF{V.LT.0.0)THEN 
SGH— 1.0 

ELSE  IF(V.EQ.O.O)TBEN 
SGN-0.0 
ELSE 

SGN-1.0 

ENDIF 

RETURN 

END 


SUBROUTINE  INTFUN(X,TIME,T,tI) 
IMPLICIT  REAL* 8  {A-H,0-Z> 

X  -  STATE  VECTOR 
TIME  -  RUNNING  TIME 
T  -  TIME  INTERVAL 
N  -  DIMENSION  OF  THE  STATE  VECTOR 
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!>IHEKSIC>i  X(20)  ,D(20;  ,A(20,5)  ,X6(20) 

DATA  AUSE/1/ 

C  INITIALIZE 

KIHC-  n 
EMIN-  0 
EMAX»  l.H-05 
HJflH-  1.2-07 
TIN-  TIME 
TOUT-  TIN  +  T 

C  TIN  -  PEGIHHING  OP  INTERVAL 

C  TOUT  -  BKD  OF  THE  INTERVAL 

IPCKUSB  .RE.  OJH-T 
IPfKOSE  ^EQ.  0)K-ESAVE 
KH-  0 
KUSE-  0 

C  IKTEGHATION  ALGORITHM  BETWEEN  22  AND  19 

DO  101  I-1,H 

101  2B(I)-  X(I) 

22  K-i 

15  CaiLi  DERIV{X,TIME,D) 

GO  TO  (100,200,300,400,SOO),K 

100  DO  11  I-i,S 

A(I,K)-  D(I)*H/3. 

11  A!I}"  2B(I>+A(I,1) 

TIME-  TIN  +  H/3.00 
X“  KM 
GO  TO  15 

200  DO  12  1«-1,H 

A(r,K>-  D(l)»H/3.00 

12  X(I)-XE(I)+0.50''{  A(I,1>+A(I,2)  ) 

K-K+1 
GO  TO  IS 

300  DO  17  1-1, K 

A(I,K)-  D(I)«H/3.00 

17  X(I)-  XB(I)+{3.0  •  A(I,1)+9.0»A(I,3))/8.00 

TIME-  TIN  +  O-SO^K 
X-  K+1 
GO  TO  15 

400  BO  la  1-1, N 

A(I.K)-  D(X)»H/3.00 

18  X(I)-XB(IJ+(3,«A(I,15-9.«A(1,3)+12.*A(I,4) )/2.0 

TIME-  TIN  +  H 
K-K+1 
GO  TO  15 

SCO  DO  19  I-1,N 

A(I,K)-  D(I)»H/3.00 

19  X(I)-XB(I)+.5*(A(I,1>+«.*A<I,4)+A(I,5)) 

C  COMPUTE  THE  TRUNCATION  ERROR 

C  INTEGRATION  ALGORITHM  BETWEEN  22  AND  19 

ERROR-  0.00 

DO  21  I-1,H 

TE-A(I,l)-{9.*A(l,3)-8.*A(I,4)+A<I,5))/2.00 
21  ERROR-  DHAX1( ERROR, DASSJTE)) 

IF(  ERROR  .GE.  EHAX  }GO  TO  33 

DO  32  I»1,H 

32  SB(I)-  X(I) 

TIN-  TIME 

if;tikb  .eq.  todt)go  to  39 
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THEM-  TOUT  -  TIME 
ip(trem  .gt.  h)go  to  31 
ESAVE>«  H 
KH-  1 
H-  TEEM 
GO  TO  22 

31  1P{TREH  .LT,  (2.00»H))GO  TO  22 

1P(EBR0R  .GT.  EHINjGO  TO  22 

KIKC-  KINC+1 

1P(KIHC  .LT.  3)GO  TO  22 

H-  2.C0*H 

KING-  0 

GO  TO  22 

33  H-  H/2.00 

IF(B  .LT.  HHIN)GO  TO  35 
TIME-  TIN 

DO  34  I-1,H 

34  X(I)-  XB(I) 

GO  TO  22 

35  HRITE(»,*) 'H  IS  LT  HMIH.  TERMINATED.' 

GO  TO  40 

39  IP(KH  .EQ.  IjRETORN 

HSAVE-  H 

RETURN 

40  STOP 
END 
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